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On Local Flat-Plate Similarity in the 
Hypersonic Boundary Layer' 


F. K. MOORE* 
Cornell Aeronautical Laboratory, Inc. 


Summary 


A study is made of Lees’ “‘local flat-plate similarity”’ rule for 
the hypersonic laminar boundary layer. It is shown that this 
rule is exact under assumptions commonly invoked in the in- 
viscid theory of hypersonic flow. 

Beginning from this theoretical basis, a modified local flat- 
plate similarity scheme is derived, involving separate rules for 
velocity and enthalpy profiles, and is compared with exact simi- 
larity solutions and with the existing theory of hypersonic lead- 
ing-edge interaction. 


Symbols 


1 + 2.60 7,,/To) 
dimensionless measure of pressure gradient (= U’;’g?) 


measure of wall temperature ( 


ratio of specific heats 

displacement thickness 

small quantity ( (gs/gs) — 1) 

similarity variable (= Y V » g) 

small constant coefficient 

inclination of surface to free stream 

constant of proportionality in pressure-viscosity 
relation 

coefficient of viscosity 

kinematic viscosity coefficient 

density 

stream function (U = y,) 

arbitrary function 

speed of sound 

arbitrary function 

constant 

skin-friction coefficient, heat transfer coefficient 

Blasius profile function 

velocity profile function 

scale function in similarity variable (¢ 


> 
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total enthalpy in boundary layer 

Mach number 

constant exponent 

coefficient in expression for overpressure due to 
strong interaction 

static pressure 

Reynolds number [Eq. (51)] 

enthalpy function 

static temperature 

transformed velocity [= ao(u/a;)|} 

velocity parallel to surface 

transformed velocity normal to surface 

M w~*®/V Rx 


coordinates along and normal to surface, trans- 


hypersonic interaction parameter ( 


formed to account for compressibility 
physical coordinates along and normal to surface 


X 
function of X ( ja U,dX ) 


be 
function of x ( f, u, pidx 


Subscripts 


stagnation conditions outside boundary layer 

conditions at outer edge of boundary layer 

evaluation at surface 

evaluation in free stream 

associated with similarity-rule modification for 
velocity profile 

associated with similarity-rule modification for 


enthalpy profile 


Superscripts 


primes denote ordinary derivatives 
. = transformation of surface coordinate under which 
local similarity applies 


Introduction 


OME YEARS AGO, Lees! put forward a simple method 
for calculation of the laminar boundary layer in 
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Fic. 1. Variation of skin friction and heat transfer rate with 
favorable pressure gradient, when U,; = CX”. 


hypersonic flow. This method, known as the local 
similarity rule, is based on the observation that the 
gas near the body surface (at the base of the boundary 
layer) tends to be much cooler and denser than the gas 
in the rest of the boundary layer. This is true pro- 
vided the surface is very much colder than the stagna- 
tion temperature of the oncoming gas. Consequently, 
it is argued, streamwise pressure gradients which may 
exist in the external flow are not very effective in 
changing boundary-layer profiles of velocity and en- 
thalpy. It then follows that the hypersonic boundary 
layer over a cold wall may be described quite well by 
the classical Blasius solution, applied to compressible 
viscous flow over a flat plate, taking account of pres- 
sure changes only insofar as local pressure affects local 
Reynolds number. This effect is accounted for by 
making a suitable transformation of streamwise co- 
ordinate. 

The local flat-plate (or Blasius) similarity rule, as 
such, does not constitute a theoretically derived result, 
but is, rather, a reasonable and very practical hypothe- 
sis which is easy to apply in a variety of hypersonic 
flow problems. The accuracy and limitations of the 
local similarity rule may be studied by making com- 
parisons with the results of calculations carried out by 
Cohen and Reshotko? in which the boundary-layer 
differential equations are solved exactly (under certain 
physical assumptions) for a class of flows subject to 
pressure gradient, having various ratios of wall tem- 
perature to stagnation temperature. Fig. 1 shows com- 
parisons of exact and similarity-rule calculations of wall 
shear and heat transfer rate as functions of pressure 
gradient. These comparisons will be made more ex- 
plicit in a later section; it suffices here to observe that 
the local similarity rule, which would give shear and 
heat-transfer rate to be independent of pressure gra- 
dient, represents heat-transfer rate fairly accurately 
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over a wide range of favorable pressure gradient if the 
wall is cold, and less well if the wall is insulated 
Representation of skin friction by local similarity is 
much less successful. 

A purpose of the present paper is to show the sense 
in which the local Blasius similarity rule follows theo- 
retically from the boundary-layer equations. In par- 
ticular, it will be shown that such local similarity is 
strictly correct under a set of assumptions commonly 
invoked in the literature of inviscid hypersonic flow 
theory, especially those pertaining to the shock layer*4 
and hypersonic small disturbance theory.®® This 
connection between inviscid and viscous theory’ is an 
important consideration in the study of hypersonic 
flows in which shock-wave and boundary-layer phe- 
nomena must be considered together. Cheng ef al.* 
have recently used this approach in the study of hyper- 
sonic flow over slender bodies with mixed effects of 
bluntness and boundary-layer displacement. 

The idea that local Blasius similarity and inviscid 
hypersonic theory are in a sense consistent may be 
made plausible in the following way. In hypersonic 
flows, inviscid streamlines tend to be very stiff, that is, 
to have such high inertia that they are not easily af- 
fected by pressure gradient. For example, in various 
modified Newtonian flow theories it is assumed that 
hypersonic streamlines are deflected only by direct im- 
pingement on a surface (and a subsequent constraint 
to follow the surface). In a similar way, hypersonic 
streamlines are laterally stiff—that is, transverse pres- 
sure gradients have small effect in causing lateral de- 
flections. It is for this reason that a tendency exists 
for streamlines to follow the geodesics of the surface. 
This aspect of inviscid hypersonic flow theory has re- 
cently been treated by Cheng.‘ 

Now, the local Blasius similarity rule assumes, in 
effect, that streamline stiffness exists within the 
boundary layer as well as beyond it. That is, stream- 
lines, after entering the boundary layer, maintain 
sufficient inertia (especially if density increases toward 
the surface) that they are not greatly affected by pres- 
sure gradient. For three-dimensional boundary-layer 
flows, local Blasius similarity implies the absence of 
secondary flow. Absence of secondary flow within 
the boundary layer is, of course, directly analogous to 
the tendency of external, inviscid streamlines to follow 
surface geodesics. The foregoing discussion is intended 
only to sketch the reasons why one may suppose that a 
connection exists between the local-similarity rule of 
Lees and the rather more well-developed theory of in- 
viscid hypersonic flow. In the analysis to follow, Lees’ 
original similarity rule will be developed in such a way 
that the comparison with Cohen and Reshotko’s work 
may be easily drawn, and so that the limiting condi- 
tions under which the local Blasius similarity rule is 
exact may be easily identified. 

A second purpose of the present paper is to propose 
a modification of the local similarity rule which takes 
into account, in an approximate way, the terms ne- 
glected in the formulation of the simple local similarity 
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rule. The modified local similarity rule will be com- 
pared with the theories of strong and weak inter- 
action!':!* for hypersonic flow over a flat plate. 
Throughout the present study, it will be assumed 
that the gas is ideal, with constant specific heats, 
that the Prandtl number is unity, and that viscosity 
varies linearly with temperature. This set of assump- 
tions is not essential for the application of local simi- 
larity rules, but is chosen for convenience, in view of 


the special purposes of this paper. 


Local Flat-Plate Similarity in Simplest Form 
Derivation of Similarity Rule 

We begin by writing the compressible laminar bound- 
ary-layer equations in the form given by Cohen and 
Reshotko,’ specializing their equations to the case of 
Prandtl number equal to unity. Under a variant of 
the Stewartson-Illingworth transformation, and an 
assumption of linear dependence of viscosity on tem- 


perature, 
, . ee 
; ap ; “ pa 1 
X= [ jd; V= [ = dy; aed - 
J0 avy J0 pollo Mo T 


the equations of continuity, momentum, and energy 
take the form 

Uy + Vy = 0 

UUxy + VUy = UUix(] + S) + wl yy (2) 

USxy + VSy = wSyy | 
where 

U = (ao/ai)u; S=(h/h) — 1 (3) 

The boundary conditions to which this system of 
equations is subject are: 


UixX, 0) = Viz, 0) = 0: UC, oe) = UX) 


Next, following Lees,' we attempt to specify that the 
solution has similarity in terms of a new variable 


C= V/V» 2(X) (5) 


where g(X) is a scale function to be determined. In 
terms of the new variable ¢, a stream function is defined 
in the same form one would use for flat-plate flow, and 
the enthalpy function also is redefined: 


vV=Vm Ui 2X) fe, X); S= SE, X) (6) 
whence, 
U= yy = Uf; 


Under these changes of variable, boundary conditions 
(4) would become 


S(O, X) = f,(0, X) = 0; fy(o, X) = 1) 
» (7) 


S(O, X) = SS; S,(@, ZX) =0 \ 


and differential Eqs. (2) would take the form 


fere + g( Ug) fee + B11 + S — f,?) 


c¢ 


l 12 ~( lex — Ixlee) (S) 
Se; a 2£\ l 4 fA Ug" ( f-Sy — fy S;) (9) 

where 
B = U'g? (10) 


Of course, if local similarity holds, f and S are func- 
tions only of ¢, and the right-hand sides of Eqs. (8) and 
(9) may be set equal to zero. 

Eq. (8) shows that, in general, a local similarity de- 
fined by f = f(¢) is impossible, because one cannot de- 
fine g in such a way that all dependence on X disap- 
pears. However, if 8 should be small for the particular 
problem under consideration, or at least if the term in 
Eq. (8) involving 8 is small, then local similarity is ac- 
complished by selecting g(X) to satisfy 


g(Uig)’ = J (11) 


Under a boundary condition which identifies x = 0 as 
the starting point of boundary-layer growth, differen- 
tial Eq. (11) may be solved to yield 


—; 
g = (1/0) y | UidX (12) 
~/ OV 


If we now return to Eq. (5), we may define the simi- 
larity variable ¢ in the form conventional for the flat- 
plate boundary layer, defining a quantity X* which 
gives the transformation of streamwise coordinate under 
which local similarity applies: 
: xX 
t= ‘ _— Y; X*= : | UX (13) 
2X . l 1 0 

Returning to physical coordinates through Eq. (1), we 
may then write 











re uy Po 7 p 
; = — f dy; x* = 
Zrx* pi JO po 
] = 
pPimdx (14) 
Pitty J 0 
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a we 
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DIMENSIONLESS NORMAL COORDINATE, ¢ 


Fic. 2. Functions appearing in similarity theory of boundary 
layer. 
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By use of Eq. (10), we find that the pressure gradient 
factor takes the form 


U,' c* 
B=2 I U\dX (15) 
l — 0 

or, returning to physical quantities, 
1 + '/e(y — 1)M;? (pum) f piuidx , 
2 V2 0 (16) 
_— (pia)? 


If local similarity is to hold, this latter quantity must 
either be negligible or, as a special case, be a constant. 
Eqs. (8) and (9) may now be written in the simple form 


B= 


+ ff" + B( + Ss = a" — 0 (17) 

5S” + fS’' = 0 (18) 

and, provided the term bearing 6 may be neglected, 

the solution may be written down in terms of the 
Blasius function 


f=F®); S=S8,[1 -— F’'O]) (19) 
where 
F'’ + FF" =0; F(0) = F'(0)=0 Fo) =1 
(20) 


Application to Falkner-Skan Flow 


We shall return in the next section to the question 
of the neglection of the pressure gradient term in Eqs. 
(8) and (17). First, however, we may note that appli- 
cation of Eq. (19) to the special case of Falkner-Skan 
flow 

U, = CX” m = (const.) (21) 
permits an assessment! of the local-similarity rule by 


comparison with the exact solutions of Eqs. (17) and 
(18) by Cohen and Reshotko? for constant values of 8. 
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Fic. 3. Comparison of skin friction for boundary layer with 
pressure gradient and similarity, as given by modified local 
Blasius similarity rule and exact calculation.? 
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Substitution of Eq. (21) into (15) yields 


B = 2m/(1 + m) (22) 
and substitution into Eq. (13) yields 
X 
im o} os r 
” . t+" “e2= 5 if (23) 
2 yyX 


These are precisely the variables in terms of which the 
solutions of Cohen and Reshotko® are expressed. 
Therefore, by adopting the local-similarity solution 
defined in Eq. (19), we are duplicating the exact solu- 
tion at 8 = 0, and assuming that solution to hold for 
a more extended range of 8. Thus, for wall shear and 
surface heat-transfer rate, respectively, local flat-plate 
similarity yields 

f’(0) = 0.470; S’(O) = 0.470[1 — (7,/To)] (24) 


We see from Fig. 1 that local flat-plate similarity is 
quite successful as regards surface heat-transfer rate 
at a cold wall, that is, 7),/7) is approximately zero. 
For skin friction, the approximation is much less satis- 
factory, and in either case the approximation suffers in 
accuracy when the wall is not extremely cold. 


Relation to Shock-Layer Theory 


We turn now to the question of when the last term 
of Eq. (17) may properly be neglected. Lees!" has 
pointed out that a cold wall tends to make this term 
small—i.e., using Eq. (19) with Sw = —1 yields 


1+S—f”? = F(1 — F’) 


This quantity vanishes both at the wall and far from 
the wall, and never has a value greater than 1/4. 
However, for the Blasius function, the remaining terms 
in Eq. (17) do not have a value exceeding 1/4 anyway. 
Thus, while a cold wall enhances the accuracy of the 
local similarity rule, the crucial quantity is @ itself. 
Eq. (16) displays a suitable expression for 8. For 
purposes of discussion, we may assume that the quan- 
tity placed in square brackets in Eq. (16) is generally 
of unit order. If the flow Mach number at the outer 
edge of the boundary layer is small, then the remaining 
factor in the expression for 8 is also of unit order. 
However, in hypersonic flow, the expression involving 
M, can be very small. If we consider that .1/,? is 
very large, 6 will be small provided 


(1/M\?) + '/x(y -— 1) <1 (25) 


This means that not only must ./,? be large, but 
'/o(y — 1) must be small. Now, various inviscid hyper- 
senic theories embody special cases of these assump- 
tions. We may mention Cheng’s treatment‘! of the 
shock-layer problem for pointed bodies, Cole’s shock- 
layer theory® for slender bodies, and a recent study* 
of slender-body flow with both bluntness and boundary- 
layer displacement. In the foregoing examples, sim- 
plifying use is made of the assumption that !/2(y — 1) 
< 1, and the accuracy of this assumption has, in certain 
cases, been verified.* + 
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Local flat-plate similarity is not specifically applicable 
near stagnation-point regions, where ./,° may not be 
large. For example, in the special case of an attached 
oblique shock, .1/,? is large only if 


M..2 cos? 6 > 1, and !/o(y — 1) tan?6<1 (26) 


and therefore the vicinity of the place where surface 
inclination (@) equals 7/2 should be avoided. Thus, 
local Blasius similarity does not form a part of shock- 
layer theory for blunt bodies,* except far downstream 
of a localized bluntness. This is not to say that local 
similarity may not be used near a stagnation point as 
has been done, for example, by Lees‘ and Vaglio- 
Laurin,® but rather that, when used in such regions, 
the validity of local similarity is not assured on theo- 
retical grounds, but must be justified on the basis of 
practical accuracy. As a further specialization of the 
requirement given in Eq. (25), we may observe that 
the value of 8 is dominated by !/2(y7 — 1) when (1 /./,”) 
is much less than '.(y — 1). This, as we may verify 
from the oblique shock relation, requires that tan’ 6 
be much less than one. In other words, flow deflec- 
tions must be small. In this case, we may make an 
explicit connection between the local flat-plate simi- 
larity rule for the hypersonic boundary layer and the 
hypersonic small-perturbation theory to which the 
assumption y — 1 <1 is added. Cole® has used such 
an inviscid theory for the study of slender shapes, 
adding the further assumption of strong shock waves. 


Extended Local Similarity Rule 


Improved Rule for Velocity Profile 

We have seen in the previous section that Lees’ local 
Blasius similarity rule for the hypersonic boundary 
layer is correct in the limit [Eq. (25)] of assumptions 
which also form the basis of a body of inviscid hyper- 
sonic theory. Since, in practical problems, '/2(y — 1) 
ranges between 10 and 20 percent, it has been necessary, 
for inviscid theory, to evaluate the accuracy of an 
assumption of y close to one by carrying out higher 
approximations in this quantity.** By analogy, then, 
to this previous body of hypersonic flow theory, it 
should be useful from the theoretical standpoint to ex- 
tend the local Blasius similarity rule to include, at 
least approximately, the effect of nonzero values of 
'4(y — 1) or, more generally, the quantity displayed 
in Eq. (25). Doing this should provide a theoretical 
measure of the accuracy of the local-similarity hy- 
pothesis in rather general circumstances. A second rea- 
son for attempting an extension of the local similarity 
rule is the more practical one of providing a rule which 
yields skin friction and heat-transfer results of greater 
accuracy than the original rule of Lees provides. The 
latter purpose requires that no significant sacrifice in 
simplicity be made in the modification of the similarity 
rule. 

We may outline the procedure which will be followed 
here in terms of Eqs. (8) and (9). We notice first 
that the pressure gradient parameter 8 appears only in 
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Fic. 4. Comparison of surface heat transfer rate for boundary 
layer with pressure gradient and similarity, as given by modified 
local Blasius similarity rule and exact calculation.* 


the momentum equation. Thus, the result shown in 
Fig. 1 is not surprising, namely, that the error incurred 
by local flat-plate similarity for the Falkner-Skan 
problem is more serious as regards velocity profile than 
enthalpy profile. Since the last term on the left of 
Eq. (8) will be considered a correction term, we are 
justified in evaluating it by use of the local Blasius 
similarity solution, Eqs. (19) and (20). We therefore 
may rewrite Eq. (8) in the form 


grrr +. [g,(Uig,)’ +. A) ff” +4- 
(B[F’(1 — F’) + (7,/T)(1 — F’)| — AFF"} = 0 


(27) 


where, again, the right side is dropped because we in- 
tend that f ~ f/(¢). The quantity 1, appearing in Eq. 
(27), is an arbitrary coefficient representing a purely 
formal addition and subtraction of the quantity //” 
within the differential equation. Observing now the 
term displayed in braces in Eq. (27), we see that it 
should be possible to choose -1 as a function of 6 and 
7\,/T) in such a way that the effect of the term in 
braces on velocity profile is suitably minimized. We 
then may determine a new scale function g, for local 
flat-plate similarity by setting the square bracketed 
expression appearing in the second term of Eq. (27) 
equal to one. Returning to the brace appearing in the 
above equation, we see that this approach should be 
particularly successful in the case of a cold wall be- 
cause the functions F’(1 — F’) and FF” have certain 
similarities. They both vanish at the surface and in 
the free stream. Fig. 2 displays these two functions 
(and another, which will appear in subsequent con- 
sideration of enthalpy similarity), calculated from the 
Blasius profile, and shows that these functions have a 
maximum difference of about ().17 within the boundary 
layer. This difference, however, has a vanishing aver- 


age across the boundary layer; that is, 


f FF'dt = f F’'(1 — F’)dé = (F’(0) = 0.470 
0 0 


(28) 


Thus, for the cold wall, if we determine 1 by a require- 
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ment that the average error represented by the curly 
brackets in Eq. (27) should vanish, we would have 
simply that 4 = £8. Including the effect of wall 


temperature, and noting that 
f (1 — F’)d¢é = 1.22 = 2.60F"(0) (29) 
0 


we find 
A = [1 + 2.60(7,,/7») |B = aB = aU ,'g (30) 


where the quantity in brackets is, for future con- 
venience, identified by a new symbol a. 

We then proceed to determine a new modified scale 
function g, by setting the term within the first square 
brackets of Eq. (27) equal to unity: 


g(Uig;))’ + alg? = 1 (31) 


For velocity profile, this equation replaces Eq. (11), 
and instead of Eq. (12), the solution 

1 — x 
f= 7124 3 U;' +a dX (32) 
gr l - + a¥ | 1 


0 
applies. We may then represent the velocity-profile 
function as the Blasius function of the new variable ¢,: 
' U; . 
f=fG): ¢= = } (33) 
2X; 
where the new coordinate transformation replacing 


that given in Eq. (13) is 


U,' ta dX (54) 


co 
x 
] 
, 


It should be noted that the new transformation [Eq. 
(34) | does not reduce to Eq. (13) for a cold wall, be- 
cause the cold-wall limit of a is one and not zero. 
Rather, it is clear that the coordinate transformations, 
Eqs. (13) and (34), can approach each other only when 
neither one differs substantially from unity. This point 
will be discussed more fully in a subsequent paragraph. 

It may be observed that since Eq. (34) is obtained 
by causing an error term to have a vanishing average 
across the boundary layer, there is a connection be- 
tween the present procedure and the Pohlhausen method 
of approximate boundary-layer analysis. One may 
say that the procedure adopted here is, in effect, a 
Pohlhausen analysis in which no profile shape changes 
are taken into account. In effect, boundary-layer 
thickness is the only parameter of the present method. 
This severe limitation is in keeping with the present 
purpose, which is simply to modify the local similarity 
rule. 

A connection may also be drawn between the meth- 
ods under present discussion and the Gortler series 
method,!*!* which is expressed in terms of the vari- 
ables appearing in Eq. (13). The leading term of 
GOrtler’s series is a “‘local-similarity’’ formula (not 
necessarily for 8 = 0, however), and succeeding terms 
result from an expansion of both y and 8 in powers of 
X*, In contrast, we consider here that 8 is small, 
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MEASURE OF DISPLACEMENT THICKNESS 


Fic. 5. Comparison of displacement thickness for similar 
boundary layer with pressure gradient, as given by modified 
local Blasius similarity rule and exact calculation,? with cold 


wall (7/7> = 0). 


but with rather arbitrary dependence on X*, the asso- 
ciated correction being obtained by an averaging proc- 
ess rather than a Taylor series expansion. 

Having modified the scale of the velocity profile 
under local similarity, we may expect a consequent 
modification of skin friction, the accuracy of which 
we may assess by comparison with the similar solutions 
of Cohen and Reshotko.? To carry out this comparison 
we recall that the original variable ¢ is exactly the 
similarity variable in which Cohen and Reshotko’s 
results are expressed [Eq. (23)]. Therefore, we first 
express ¢,; in terms of ¢ by substituting Eq. (21) into 
Eq. (34), and then comparing Eqs. (33) and (23). We 


find 
cf, = ¢V1+ of (35) 


showing that, for this problem, the modified velocity 
profile has a steepness which increases with 8 according 
to the factor V1 + a8. Using this result, we find a 
formula for skin friction to replace that shown in Eq. 
(24): 


Uf") le = F"(0)V1 + a6 = 
0.470V1 + [1 + 2.60(T,/T>) 18 (36) 


Fig. 3 shows the application of Eq. (36) to the Falkner- 
Skan problem, and a comparison with Cohen and 
Reshotko’s exact results. It is seen that the agree- 
ment for a cold wall is very good, and, for an insulated 
wall, is quite satisfactory in terms of the correction re- 
quired for complete agreement. We may note that 
the modified rule is quite accurate for plane (6 = 1) 
and axisymmetric (6 = 1/2) stagnation points, as well 
as for moderate adverse pressure gradients. 
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Improved Rule for Enthalpy Profile 

Fig. 1 indicates that the original local-similarity rule 
[Eq. (14)] gives surface heat transfer fairly accurately 
for an extended range of 8; not as accurately, however, 
as skin friction is given by the improved rule. For this 
reason, and also because the determination of displace- 
ment thickness will appear to be sensitive to the ac- 
curacy of the rule for enthalpy, we turn now to the 
approximate solution of Eq. (9) for S. 

First, we observe that the new transformation X ,* 
[Eq. (34)] is not appropriate for enthalpy. The rule 
for enthalpy should be much nearer the original rule 
[Eq. (13)]. This may be inferred either from Fig. 1, 
or from Eq. (9). We see that the departure from the 
simple rule obtained by setting g(Uig)’ = 1 is due only 
to the fact that the coefficient f is a function of a vari- 
able ¢, which is presumably slightly different from the 
variable, which we shall call ¢,, appropriate for S.  Ex- 
pressing /(¢,) as a function of ¢, and X, we find, from 
Eq. (33), to a first approximation in (g,/g,) — 1, an 


expression 
f F(t.) + [(ge/gs) — 1)f.F’(S.) — F(E)| +... 


which we may substitute into Eq. (9), noting that since 
f(€., X), the term fyS; on the right of Eq. (9) must 


$ 
be retained as a small correction term. In the various 


terms involving the small quantity « = (g./g,;) — 1, 
we replace S’ by its approximate value, —S,,F"(¢). 


Also, just as in Eq. (27), we add and subtract a term 
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Fic. 6. Comparison of displacement thickness for similar 
boundary layer with pressure gradient, as given by modified local 
Blasius similarity rule and exact calculation,? with insulated 
wall (74/79 = 1). 
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Bf S;, finally obtaining the following equation: 
S” + [g.(Uig.)’ + B|FS’ + 
(—S,,) j le + Uig.2e’ |¢F' PF” — 
[e+ Uig.2e’ + BJFF’; = 0 (37) 
As in the analysis following Eq. (27), we may choose 
B so that the average of the term in braces vanishes. 
The function ¢F’F” is similar in shape to FF”, as may 
be seen from Fig. 2, and has an integral across the 
boundary layer equal to 1.80 F”(0). Thus, we set 
B 0.80(e + Uig.?e’) (38) 
and establish g, as the solution of the equation [analo- 
gous to Eq. (31) | 
g.(Uig,)’ + O.80(e + Uig,?e’) l (39) 
Eq. (39) may be solved, under the assumption that the 
difference between g, and g is small, to yield 


fi -4(4 | 
g, = g}! — * ~ii> 
0. sf" g , : 
— | —1)2Z-2/9dZ (40) 
Sl 7 VU Ly 


“2 [ Uidx (41) 
J 0 


CX", Eq. (40) gives 


where 


For the Falkner-Skan case, LU 
the result analogous to Eq. (35): 


o% = ¢[1 + (2/7)(V1 + a6 — 1)] (42) 


The quantity in square brackets gives the factor 
whereby surface heat transfer increases with 8. The 
comparison of the increase given by the modified en- 
thalpy rule [Eq. (42)] with the exact calculations is 
shown in Fig. 4. The agreement is quite satisfactory, 
especially at the axisymmetric stagnation point (8 
1/2). 

We may note that Eq. (39) may easily be solved 
Results show that 
1), the 


exactly in the Falkner-Skan case. 
in the most extreme condition (8 = 1, 7/7» 
correction term (2/7)(W1 + a8 — 1) in Eq. (42) is 
in error by only about 15 percent; this tends to justify 
the accuracy of Eq. (40) as a solution for g,. 


Formulas for Boundary-Layer Quantities 
Returning now to physical coordinates, we may write 
the coordinate transformation for the modified Blasius 
similarity rule for velocity [Eq. (34) ] in the form 
] 


ex 
, | pit M,** dx (43) 
/J 0 


Xy = 9 
pit. M2" 


which, when the velocity profile is under consideration, 
should replace the quantity x«* appearing in Eq. (14). 
Again, we see that there is an essential difference be- 
tween Eq. (43) and Eq. (14) which does not disappear 
in the limit of a cold wall. However, it may be shown 
from inviscid hypersonic theory that the Mach number 
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at the outer edge of the boundary layer would tend to 
vary with x much more slowly than pressure. In fact, 
for particle-isentropic flow at the outer edge of the 
boundary layer, the variation of Mach number would 
vanish in comparison to pressure variation, in the limit 
of the assumption given in Eq. (25). 

Similarly, for the modified enthalpy rule, Eq. (40) 
yields 

! 


4 x* 
x.* ae x* j— — | 4. 
. | aes ) 
10 - 2(x) x* 9/0 
a? e eg g "ds (44) 
0 Xy 


= f uypidx (45) 
0 


We now proceed to write down the appropriate 
formulas for certain boundary layer quantities of in- 
terest. Under the modified local flat-plate similarity 
rule, the heat transfer coefficient is given by 


where 


De] 


u (SY) x Av 2 . 
gw =* = 9.470 — _-* we 
pil) (—S,,) Pi 2x, * 1 Po 
while, under modified flat-plate similarity, the local 
skin friction is given by 
w\Uy)w 2ry y _ 
Hoo My - = 0.470-— - Pi (47) 
(1/2) pitts” pl UX ,* po 


Cy = 


The displacement thickness involves both the ve- 
locity profile and the enthalpy profile, and its evalua- 
tion involves an application of both rules. Beginning 
with the basic definition of displacement thickness, we 
find, after some manipulation, the expression 


t= | [1 — (pu/ pity) |dy 
0 


0 2r (" ing / — | 
pi Vox pi ‘ —" 4 9 384 Y M2 + 
Po 7\ 2 


= 1,22 


pi ut} 


(WV (x,*/x,*) — 1][1 + 1.38 2 = 


ult (48) 


In the special case of small deflections in the outer 
flow, for which '/o(y — 1) M,? > 1, Eq. (48) may be 
rewritten: 

po . /2rAx* pi y — 1 


6* = M? xX 
ey uy Po 2 


iia ee ag ar 

oe 1 + 2.60 7 a = g.00 
2° x,*\ ]) 1% 
xt y* { (4) 


This expression may be compared, for the further spe- 
cialization l) = CX”, with the exact calculations of 
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Cohen and Reshotko.? Figs. 5 and 6 show the 


brace of Eq. (49), for a cold wall and an insulated 
wall, respectively, and the corresponding factor calcu- 


lated from reference 2. Under the original local simi- 


larity rule,’ this quantity would be constant. The 


modified similarity procedure is seen to be quite accurate 
for the insulated wall, but much less accurate for the 
cold wall. This result, at first sight surprising, may be 
inferred from Eq. (49); when 7;,/ 7% (), the difference 
of scales for enthalpy and velocity profiles enters with 
a factor (3.6) quite large compared with unity. The 
result is then very sensitive to the accuracy of estima- 
tion of scales, and to the neglection of changes in pro- 
file shape. When 7),/7, = 1, on the other hand, the 
result is not so sensitive to scale estimation. 
Comparison with Hypersonic Flat-Plate Theory 

We have, in the previous section, and in Figs. 3-6, 
already assessed the accuracy of the Blasius similarity 
rule when applied to Falkner-Skan flow. We may 
provide an interesting further test of local similarity 
theory by comparison of the present results with the 
well-developed theory of viscous hypersonic flat-plate 
flow. In this problem, the effect of boundary-layer 
displacement in raising the pressure on the plate and 
modifying skin friction and heat transfer is assessed. 
It has been necessary to consider two extremes of this 
problem. One, the so-called weak interaction case,!! 
considers the boundary layer and associated inviscid 
flow on the flat plate far downstream from the leading 
edge, where the action of the displacement thickness 
is quite weak. Nearer the leading edge of the flat 
plate, strong interaction”: | holds. 


Weak Interaction 

We first apply the local similarity theory to the case 
of weak interaction, and then compare with the exact 
perturbation solutions given in reference 11. An 
iteration procedure is appropriate here, in which, to a 
first approximation, the external velocity and pressure 
are constant along the plate. To this approximation, 
the displacement thickness takes its classical form 


y—-1M = ( T 
6* = 0.664 | 2.60 *) (50) 
v4 VR, + T 5 


which may, of course, be written down directly from 
Eq. (49), noting that conditions at the outer edge of the 
boundary layer are effectively the same as those in the 
free stream, making the hypersonic small-perturbation 
approximation that 1/, is equal to J/., setting all 
x*’s equal to x, and making the identification 


R, = (mx/Av)(T0p~ / pol =) (51) 


If the displacement thickness interacts only weakly 
with the outer flow, one may write 


pi/p> = 1 + yM,.(06*/Ox) +... (52 


or, upon substitution of Eq. (50) 
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D1 (7 — 1) Pa 
P - 1+ 0.664 (ci + 2.60 — “ie re 
pb } To 
(53) 
where 
x = V.°/V R, (54) 
The problem now is to assess the effect of the changed 
pressure level given by Eq. (53) on skin friction and 
surface heat transfer. To do this according to local 
flat-plate similarity, one must first evaluate the trans- 
formations to x,*, x,*. In inviscid hypersonic small- 
disturbance theory, it may be shown that the outer 
velocity tends to remain constant, so that we may 


write, in place of Eqs. (14) and (43), 


1 eX T,® ™ x by 
x* x | pdx; xf ~~ | PL ae (58) 
Pi 0 pi o Ti 


To carry out the indicated integrations, we observe 
that pressure and temperature have the following form 
Pp» =~1l+tn/Vxt+... / 


= 1+ [(y — 1)/ye/Vxt---) 


(56) 
7, 7 
where, for convenience, we temporarily represent the 
x-dependence implicit in Eq. (53) by a constant co- 
efficient » which is small, in the present instance of 
weak interaction. The temperature profile has a cor- 
responding form determined by the condition of isen- 
tropic flow along streamlines at the outer edge of the 
boundary layer. Inserting Eqs. (56) into Eqs. (55) 
and (44), we find 
x* =x* = x(1 +9/Vxt+...) / 
(57) 
x =x}1+ [1 — aly — 1)/yln/Vxt+...}) 
showing that, in the present problem which falls within 
the scheme of hypersonic small perturbation theory, 
the modified local similarity differs from the original 
local similarity in proportion to the quantity 1 — 
a(y — 1)/y¥. 
For use in the formula for surface heat-transfer rate 
given by Eq. (46), we write 
( -* — (]/x)f1 + O(n?)] 58 
(Pi/ Po) /X (1/2)41 + Ue") I (93) 
from which we derive 
/ 
C, = (0.332/V R,) [1 + 0(x?)] (59) 
Thus, the original and modified similarity rules both 
give the correct result!! that heat transfer rate is un- 
affected to first order in the interaction parameter 
x. According to weak interaction theory, however, 
skin friction is affected to first order. To determine 
this effect by the present modified local similarity 
scheme, we first write 
Pi/ po ] : Heat ' 
— =—|l+ae n/Vx... (60) 
Xz Xx Y 
which when inserted in Eq. (47) yields, after some 


manipulation, 


T 
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P 0.664 | y- ( T.\ if? 
= ‘ 0.332 + 2.60 — 
ver’. st 2.00 7°) - 


Kt. (61) 


It is interesting that Eq. (61) is precisely the result ob- 
tained by the solution of the appropriate exact dis- 
turbance equations for the boundary layer.'! The 
original local similarity rule would set the correction 
term proportional to x in Eq. (61) equal to zero. The 
problem of weak interaction in hypersonic flat-plate 
flow is certainly a problem to which local similarity 
should apply very well, because the pressure gradient, 
and hence 6, is assumed at the outset to be small. 
Indeed, we find that the modified local similarity rule 
has a balanced accuracy for this problem in the sense 
that, to first order in x, both heat transfer and skin 
friction are given correctly. 
Strong Interaction 

We turn now to the problem of hypersonic flat-plate 
flow near the leading edge, where the shock wave 
generated by the displacement effect dominates the 
boundary layer. In the strong interaction problem, 
Li and Nagamatsu!” showed that a similarity law holds. 
The similarity is of the form studied by Cohen and 
Reshotko,? and solutions can be constructed for the 
strong-interaction problem by use of their results. 
However, in the strong interaction problem, a severe 
pressure gradient is produced by the displacement 
thickness. Therefore, in carrying out a 
interaction analysis by means of local flat-plate similar- 


strong- 


ity rules, one is, in effect, subject to the limitations of 
accuracy cited in Figs. 5 and 6. 

Following Lees,'” we employ the 
approximation in order to relate pressure distribution 
along the plate to displacement thickness*: 


tangent-wedge 


(pi/po) = [y(y + 1)/2)M..2(db*/dx)? (62) 


It may be readily verified, from Eqs. (49) and (62), 
that consistency requires 

6* x 43/4 (63) 
Making use of Eq. (63), we can substitute Eq. (49) 
into Eq. (62) to yield an expression for pressure dis- 
First, we must determine the pertinent 
We find 


tribution. 
transformations from Eqs. (44) and (55). 


w/e = 2 (64) 


x 
x,* ] -(1-«7—"*) q j~2 7 j 
= xpi , Pi Y ax = 
x* 2 f J0 


a: | 1 
(; +a : ) (65) 
Y 


x,*/x* = 1+ (18/35)[V1 + a(y — 1/7) — 1] (66) 


which we may use in the following formula derived 
from Eqs. (49), (51), (54), and (62): 


* We do this to facilitate comparison with his results. As 
previously noted, for theoretical consistency, we should use the 
shock-layer theory of Cole® or Cheng® for this purpose. 
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3 . 
: aie (0.664) (y — 1)V vy -- 1) xX 
Pes < 


(1 + 2.60 =) x 3.60 ( es e) 
= 7) V. ‘i Ve \ ag ° 


(67) 


One may, of course, proceed to apply Eq. (67) to calcu- 
late heat transfer and skin friction from Eqs. (46) and 
(47). 

We now compare Eq. (67) with results obtained under 
the original local-similarity rule!’ and by exact calcu- 
lation.” In effect, we have only to recast the com- 
parisons shown in Figs. 5 and 6 for the strong-interaction 
problem. We write the pressure formula in the form 


pi/po = Ply, (Tw/To) Ix (68) 


and evaluate the quantity P by the three means out- 
lined above. We may note that 6 for this problem is 
simply (y — 1)/y, which is, for air, 0.286, and for 
helium, 0.4. Table 1 shows values of P calculated for 
these two gases. It is seen that for a cold wall, the 
modified local-similarity scheme underestimates the 


TABLE 1 
Values of P Calculated for Air and Helium 


Modified Local 





Original Local 








Tw/To Similarity!” Exact? Similarity 
Air He Air He Air He 
0 0.18 0.30 0.14 0.21 0.12 0.18 


l 0.66 1.10 0.47 - 0.46 0.69 


interaction effect by about 15 percent for both air and 
helium, while the original similarity rule overestimates 
by 29 percent for air, and 43 percent for helium. For 
an insulated wall, the modified rule is much more ac- 
curate than the original rule, at least for air. These 
errors are in proportion to those presented earlier in 


Figs. 5 and 6. 


Conclusions 


A study has been made of the local flat-plate simi- 
larity rule put forward by Lees! for the study of hyper- 
sonic laminar boundary layers. It is shown that the 
rule in original form [Eq. (14)] becomes theoretically 
valid when the Mach number at the outer edge of the 
boundary layer becomes large and the specific-heat 
ratio approaches unity. Thus, the similarity rule is 
the viscous counterpart of inviscid hypersonic theory 
for attached shock layers,* including hypersonic small- 
perturbation theory in the form studied by Cole.* In 
these inviscid theories, the usefulness of the assump- 
tion of specific-heat ratio near 1 has been made specific 
by calculation**® of corrections of first order in y — 1. 
Accordingly, in the present study, parallel improve- 
ments of the local-similarity rule for viscous flow have 
been sought. 

The desired correction to the original rule is obtained 


SCIENCES OCTOBER 1961 


by an averaging process, yielding separate rules for 
velocity profile [Eq. (43)] and enthalpy profile [Eq. 
(44)]. This modified scheme is evaluated in special 
cases by comparison with the exact calculations of 
Cohen and Reshotko,? and appears to be especially 
accurate, for a cold wall, in estimating skin friction 
and heat transfer, even at stagnation points and in 
moderate adverse pressure gradients. Displacement 
thickness is particularly sensitive to the essential de- 
fect of any flat-plate similarity rule; namely, that the 
effect of pressure gradient on profile shape is neglected 
It turns out that the modified local-similarity rule 
gives displacement thickness more accurately for an 
insulated wall than for a cold wall. 

The present modified local-similarity scheme is ap- 
plied to the hypersonic problem of leading-edge shock 
and boundary-layer interaction. For weak inter- 
action, the exact first-order formulas for skin friction 
and heat transfer are exactly reproduced. For strong 
interaction, modified similarity gives the over-pressure 
with an accuracy of about 15 percent for a cold plate, 
and more accurately for an insulated plate. 
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The Effectiveness of Radiating Fins 
With Mutual Irradiation 


E. M. SPARROW,* E. R. G. ECKERT,* ann T. F. IRVINE, JR.** 
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Summary 


Finned surfaces are widely used for heat exchangers which 
transfer heat by convection. This fact suggests an investigation 
of the advantages of fins for heat exchangers transferring heat 
by radiation. The analysis of this situation becomes involved 
when the fins are arranged in such a way that they “‘see’’ each 
other so that they interchange radiation. The present analysis 
formulates this problem as a system of integro-differential equa- 
tions for gray and diffusely reflecting surfaces. Solutions have 
been obtained by a weighted iteration procedure on an elec- 
tronic computer for a heat exchanger element consisting of a tube 
with external, longitudinal fins equally spaced at angles of 45, 
For the emissivity of the fin surfaces, 
It is assumed that 


60, 90, or 120 degrees. 
the values 1.0, 0.75, and 0.5 were selected. 
the finned tube radiates freely into a space at zero temperature 
and that the tube radius is small relative to the height of the 
fins. Additionally, the fins are assumed to be very long in the 
longitudinal direction. The results can also be applied to any 
system of plane radiating fins, any pair of which share a common 
edge. Temperature distributions, local and overall heat loss, 
and fin effectiveness are presented and the conditions for which 
the weight of the fins becomes a minimum are determined. It is 
demonstrated that a parameter similar in form and value to one 
used for convective fins determines optimum conditions. A few 
comparisons of some specific fin arrangements conclude the paper. 


Symbols 

Ap = profile area of a fin 

B(x) . , : 

Bly) = combined radiant flux (emitted plus reflected ) leav- 

‘ ing a position x or y per unit time and area 
(radiosity ) 

F = angle factor 

H = incident flux per unit time and area 

JC = dimensionless incident flux, 1/a7)' 

k = thermal conductivity 

i = fin height 

l = dimension normal to plane of Fig. I(a) 

N, = conduction parameter: L?o7)°/kt, when heat is 
transferred from both fin faces; L?a07)*°/2kl, when 
heat is transferred from one face 

O = overall rate of heat transfer from one face of fin 
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Or rate of heat transfer from x = O tox 
q = local heat transfer rate per unit area 
5 = surface area 

i = absolute temperature 


t = half thickness of fin 
coordinates measuring distance along fin from bas 


£9 = 
surface 

a = absorptivity 

B = dimensionless radiosity, B/oT)' 

> = opening angle 

n = fin effectiveness, Q/o7,'L sin (7/2 

4 dimensionless temperature, 7/7» 

emissivity 

r = dimensionless coordinate, y/L 

E = dimensionless coordinate, x// 

p reflectivity 

a Boltzmann's constant 

¢ = angle defined in Fig. I(a) 
Subscripts 

b = base surface 

opt optimum value 


Introduction 


aor TO SPACE TECHNOLOGY have generated 
renewed interest in radiation heat transfer. Con- 
vection is absent in atmosphere-free space and radiation 
is usually the only means by which waste heat from 
power plants, electronic equipment, or other sources in 
a space vehicle can be disposed. The surface areas 
required for this purpose are often huge, since the 
Stefan-Boltzmann law sets an absolute upper limit to 
the heat that can be transferred per unit time and sur- 
face area. Optimization of heat exchangers for radi- 
ative heat transfer is, therefore, a very contemporary 
problem. Other space-oriented applications involve 
the capture of solar radiation and the conversion of 
this energy into more usable forms. 

In examining techniques by which radiation heat 
transfer to or from bodies may be improved, it is nat- 
ural to consider the use of extended surfaces—.e., fins. 
The use of extended surfaces to augment convective 
heat transfer is common practice, and the character- 
istics of such surfaces have been studied quite thor- 
oughly. In contrast, only a modest beginning has been 
made toward analyzing the heat-transfer performance 
of radiating fins. Thus far, primary consideration has 
been confined to single radiating fins which exchange 
heat with space but do not interact with neighboring 
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(0) ANALYZED FIN CONFIGURATION 








A 





(c) FINS ON A PLANE SURFACE 


(b) FINNED CYLINDER 


Fic. 1. Schematic diagram of fin configurations. 


fins.'~* The qualitative behavior of such simple fins 
is similar to that of convection fins in that the greatest 
heat transfer takes place near the fin base (region of 
highest temperature), with less heat transfer occurring 
near the tip. The interaction of a single fin with its 
base surface has been mentioned in reference 5, and an 
incomplete analysis of the problem is reported in refer- 
ence 8.* 

The present investigation concerns the situation in 
which the radiating fins can exchange heat with each 
other as well as with the external environment. Here, 
the problem becomes more complicated with regard to 
both mathematical formulation and physical interpre- 
tation. The mathematical representation of the com- 
plex process of radiant interaction and heat conduction 
requires an integro-differential equation’ rather than a 
differential equation such as arises for the single radi- 
ating fin. By physical reasoning, it is possible to arrive 
at certain interesting conclusions about the effects of 
the radiant interaction between the fins. 

To visualize the problem, consider Fig. 1(c), which 
shows a plane isothermal surface upon which corrugated 
longitudinal fins have been erected. We want to com- 
pare the heat loss from the unfinned plane surface with 
the heat loss which will occur when the fins are used. 
First, suppose that the fins and the base surface are 
both black-body radiators (« = 1). If the thermal con- 
ductivity of the fin were infinite, then the temperature 
throughout the fin material would be identical with 
that of the base surface. Under these conditions the 
heat loss from the fin surfaces would be identical with 
that from the unfinned base surface (because the dashed 
envelope area is equal in size to the base area). If the 

* Added in proof: the fin-base surface interaction has recently 
been analyzed in reference 12. 
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surfaces were maintained black but the fin conductivity 
were finite, then the heat loss from the fins would 
actually be less than that from the unfinned base sur 
face. On the other hand, for nonblack surfaces (€ < 1), 
fins of infinite thermal conductivity would increase the 
heat loss; but this increase would diminish and would 
finally be wiped out as the conductivity of the fin mate 
rial decreases. Therefore, the use of longitudinal fins 
on a plane surface is most advantageous in situations 
where the emissivity is low and the thermal conduc 
tivity high, but the quantitative limits remain to be 
determined by analysis. 

If we shift our attention to Fig. 1(b), which depicts 
longitudinal fins on a cylindrical surface, we find that 
somewhat different conclusions apply. Here, for black 
surfaces and infinite fin conductivity, the use of fins 
will increase the heat loss by the ratio of the area of the 
dashed envelope surface to the area of the base sur 
face. For finite thermal conductivity and black sur 
faces, fins may still be advantageous, but to a lesser 
extent with decreasing thermal conductivity. For 
sufficiently low thermal conductivity, the fins may 
diminish the heat loss relative to the unfinned base sur 
face. When the surfaces are nonblack (e€ < 1), high- 
conductivity fins will definitely increase the heat loss; 
but again, this gain diminishes with decreasing con 
ductivity. Therefore, for the configuration of Fig. 
l(b), fins may be advantageous for surfaces of both 
high and low emissivity, provided that the thermal 
conductivity is sufficiently high. Thus, the use of 
fins appears to be potentially more advantageous for 
the configuration of Fig. 1(b) than for that of Fig. 1(c). 
We can generalize this conclusion by noting that the 
use of fins for increasing the heat loss becomes more at- 
tractive with increasing values of the ratio of the area 
of the envelope enclosing the fins to the base surface. 

With these remarks as general guideposts as to the 
utility of radiating fin ensembles, we turn to the quanti- 
tative determination of fin performance. The par- 
ticular configuration chosen for study is shown in Fig. 
l(a). Here are pictured a pair of longitudinal fins of 
rectangular profile which extend to a relatively great 
length in the direction normal to the page. The fin 
surfaces are gray, diffuse radiators. The configuration 
was originally selected with the longitudinally-finned 
cylinder, Fig. 1(b), in mind. To reduce the already 
large number of parameters, consideration has been 
given here to exchangers with fin height 1 much greater 
than the cylinder radius. It is clear that the results 
obtained for the system of Fig. 1(a) will also apply to 
corrugated fins on a plane wall, as shown in Fig. 1(c). 
The analysis is carried out for the condition that the 
radiation impinging on the fins from the external en- 
vironment is negligibly small, a situation which holds 
for sufficiently high surface temperatures. The specific 
role of external radiation depends on particular condi- 
tions such as the angle of orientation with respect to 
the sun and the absence or presence of focusing mirrors. 
These conditions can be included when specific designs 


are being studied. 
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EFFECTIVENESS OF 


In an earlier note® a procedure for setting up the gov- 
erning equations for radiating-conducting fins with 
mutual irradiation was outlined. Readers may find 
this prior paper a useful corollary reference source. 


Analysis 


Energy Conservation and Radiant Flux Balance 

The starting point of our analysis of the radiating- 
conducting fin ensemble is the basic law of energy con- 
servation. The application of this principle is facili- 
tated by reference to Fig. l(a), where attention will be 
directed to the energy transfers to a typical cross- 
hatched element. This element corresponds to the 
situation where there is a symmetrical heat transfer 
from both surfaces of the fin—e.g., Fig. 1(b). If one of 
the faces of the fin is essentially insulated—e.g., Fig. 
1(c)—then the element would span the entire thickness 
of the fin. Under steady-state conditions, energy con- 


servation requires that 


(dQrad net — (dQ cond.)net = OU (1) 
The net conduction is found by application of Fourier’s 
Law*: 
d d1(x) d?7(x) 
(dQcond.)net = — kil dx = —kt —— £5. 
dx dx dx? 
(2) 


where / represents the dimension normal to the plane of 
the figure and dS, = / dx is the surface area of the ele- 
ment. Next, to characterize the net radiation we 
introduce two new quantities. The first of these is the 
radiosity B, which represents the radiant flux per unit 
area and time leaving the surface; the second is H, 
which denotes the radiation incident on unit surface 
area. In general, both B and H/ will vary with posi- 
tion (x) along the fin. With these, the net radiation 


may be written as 
(dead Q)net = [B(x) — H(x)] dS, (3) 
By use of (2) and (3) the energy balance (1) may be 
evaluated, giving 
@?T(x)/dx? = (1/kt) (B(x) — H(x)] (4) 


As it stands, Eq. (4) contains three unknowns: the 
temperature 7° and the radiation fluxes B and //; and 
so, other relations between these quantities must be 
found. 

The first of these is a balance of radiant energy, which 
derives from the fact that the radiation B leaving the 
surface is the sum of the direct emission of the surface 
plus the reflection of the incident energy. The energy 
emitted per unit area is given by the well-established 
relation eo7*; while the reflected incidence is simply 
pil, where p(= | e), is the gray-body reflectivity. 


*In writing the net conduction, it is assumed that the fin 
length Z is much greater than its thickness, and therefore T is 
considered a function only of x. This approach is standard in 


fin theory. 
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And so, we can write 
B(x) = eo74*(x) + pll(x) ») 
Next, it may be observed that the energy // incident 
on @S, must be related to the radiant flux which leaves 
positions on the opposite fin. Consider atypical surface 
dS,. The radiant energy leaving dS, in all directions 
is 
B(y) dS, ba) 
Of this, an amount 
B(y) dS, dF,_, (Ob) 
arrives at dS,, where dF,_, is the angle factor of dS, as 
seen from dS,. But, from the reciprocity between 
angle factors, 
B55 BF pie = BS, Oh gin 6c) 
The definition and evaluation of angle factors is given 
in textbooks on radiant heat transfer (see, for instance, 
reference 11, p. 396). With Eq. 6(c), expression (6b) 
becomes 


B(y) | ae (Od) 


Hence, the energy incident per unit area at dS, due to 
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radiation leaving dS, is 


B(y) dF,_, (6e) 


But, dS, receives energy from all positions on the op- 
posite fin. The total contribution is obtained by in- 
tegrating Eq. (6e), and is //(x): 


y=L 
H(x) = f B(y) dF,_, (6f) 
y=0 
It still remains to evaluate the angle factor /,_,. For 
the situation under consideration here, where the fins 
extend to essentially infinite length in the direction 
normal to the plane of Fig. 1, the angle factor is given 
by the remarkably simple relation (reference 10, Eq. 
31-58) 


dF,_, = '/. d(sin ¢) (7a) 


where ¢ is the angle between the normal to dS, and the 
line connecting dS, and dS,. From Fig. 1(a), it is 
easily seen that 
: x— ycos 7 . 
sin g = a a ——— (cm) 
[(x — y cos y)? + (ysin y)?] 
With Eqs. (7), the incident flux 7 as expressed by Eq. 
(6f) becomes 


1 & xy(1 — cos? vy) 
H(x) = f ot a . — 7,dy (8) 
ost [x? + y? — 2xy cos y]”* ~ 


The energy conservation equation (4) together with 
the radiant flux balance (5) and (8) provide a complete 
mathematical description of the radiation-conduction 
process for the fin ensemble of Fig. l(a). To determine 
the three unknowns, there are now three equations. 
To complete the statement of the problem, it remains 
to state the boundary conditions, which are prescribed 


as follows: 


x = 0 (fin base) a, = De 


ae? ag (9) 
x =L (fin tip) dT/dx = 0 


The second boundary condition is an approximation 
(standard for convective fins) of the energy balance at 
the fin tip,’ and embodies the assumption that the heat 
loss from the tip is negligible compared to that from 
the fin face. 

It is useful to make the governing equations dimen- 
sionless and thereby exhibit the independent param- 
eters. Utilizing the following definitions 
8 = B/ol,', K = H/ol,', 0 = T/T, 

£=2/L,\ = 9/L (10a) 


N. = (L2/kt)oT,? (10b) 


Eqs. (4), (5), (8), and (9) may be rephrased as 


BE) = *(E) + (1 — ©)5C(E) (11) 
d*0(&)/dé* = N,[B(—) — H(E)] (12) 
en |} yy —frCl = cos? y) dd ‘ 
ete! Mak Foe oor ee 
6(0) = 1, d6/dé(1) = O (14) 
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Inspection of this set reveals that there are three inde- 
pendent parameters: the opening angle y, the emissivity 
e, and the group N,. This latter parameter charac- 
terizes the role of the heat conduction in establishing a 
temperature gradient from the fin base to the fin tip. 
It will be designated here as the conduction parameter. 
Clearly, for infinite conductivity VN, = 0, and the 
temperature is uniform over the entire fin. With finite 
conductivity, V, will be greater than zero and the fin 
temperature no longer uniform. 
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Because of the symmetry of the situation, it is easy 
to see that the variation of 6 (or B) is the same along 
both of the fins. So, @(g) and B(A) are the same func- 
tions; only the independent variables have been 
changed. Further inspection of Eqs. (11), (12), and 
(13) shows that the unknown § appears under the in- 
tegral sign as well as in other parts of the equation. 
This fact, coupled with the appearance of the deriva- 
tive d°6/dé* in Eq. (12), signifies that we are dealing 
with integro-differential equations. Since the three 
variables are intermingled in the governing equations, 
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EFFECTIVENESS 


it is evident that simultaneous solution is required. 
It might appear that the mathematical complexity of 
the problem would be reduced if 3C were eliminated 
from Eqs. (11) and (12) by use of Eg. (13). This 
would very likely be true if an analytical solution were 
possible. However, since it was necessary to obtain 
solutions by numerical means, there would have been 
no essential gain in combining the equations. 
Simultaneous solution of the governing equations 
(11), (12), and (13) yields the temperature distribution 
@ = 77/7) and radiation fluxes B and #7 as functions of 
position on the fin. As will be shown later, a knowl- 
edge of B and /Z immediately provides the desired 
information on the heat transfer performance of the 


fin. 


Solutions 


The procedure used to solve Eqs. (11), (12), and (15) 
was essentially one of iteration. The actual numerical 
calculations were carried out on a Remington Rand 
1103A (Univac Scientific) electronic digital computer. 
For a particular case the values of the parameters 7, €, 
and NV, were fixed, and trial values of 6 covering the 
range ( < \ < 1 were selected. Then, for a fixed value 
of &, the integration on the right of Eq. (13) was carried 
out numerically; this provided 3((£). The integration 
was repeated for all the mesh points in the range 0 < 
&< 1. By use of these 3C(€) in conjunction with the 
original set of 8 values, the differential equation (12) 
was solved numerically for 6(€). Then, with @(&) and 
5C(€) as input, a new set of 8 values was calculated from 
Eq. (11). To continue the iteration, these new 6 
could be used either directly as input to Eq. (13), or 
else, to avoid instabilities, they were averaged with the 
earlier ones and then inserted into Eq. (13). In either 
event, the process just described was repeated until 
convergence was achieved. 

Solutions were obtained for opening angles y of 45°, 
60°, 90°, and 120° for e values of 1.0, 0.75, and 0.5. 
For each y, « combination, the parameter ', was as- 
signed values between (0 and 4. In all, 156 cases were 
considered. It may be noted that all of the opening 
angles selected are integral divisors of 360°, and hence 
the results are applicable to longitudinal fins uniformly 
arranged around a cylinder. Emissivity values lower 
than 0.5 were not considered because it was thought 
unlikely that low emissivity surfaces would be used when 
the aim is to increase the heat loss. 

Numerically speaking, the evaluation of 3C(£) from 
Eq. (13) becomes uncertain as £ and 9 approach zero, 
since the integrand is indeterminate. It is possible to 
bypass this difficulty by physical arguments, from which 
it may be shown that 


— sin (y — oa 


2 


1 
H(0) = B(O) | 


where the bracket is the angle factor of the opposite 
fin as viewed from x = 0. Using this in conjunction 
with Eq. (11) and noting that @ = 1 at x = 0, it follows 
that 


OF 
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Consequently, the 6 and H values at x = é 0 are 
known a priori and need not be included in the iteration 


procedure. 


Heat Flux Equations 

The major result of engineering interest is the heat- 
transfer performance of the fin ensemble. The net 
local heat flux q leaving a unit area of fin surface is 


g=B-dH (15) 
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and the overall heat loss (per unit width) Q from one 


face of the fin follows as 


L O 1 
QV = f (B — H)dx or —— = f (8—35) dé (16) 
0 Cer." 0 


Since 6 and 3 are functions of y, «, and NV,, so also is 
Q/Lol;,'. In keeping with the practice followed for 
convective fins, it may be useful to rephrase the heat 
transfer results in terms of an effectiveness n, which is the 
ratio of the actual heat loss from the fin to the greatest 
possible heat loss. Now, for radiating-conducting 
fins, the greatest heat loss will occur under the condi- 
tions of a black surface (« = 1) and infinite thermal 
conductivity (V, = 0). For these conditions, the net 
loss from the fins is equivalent to that of a black surface 
stretched tightly between the tips of the fins. Conse- 
quently, for one face of a fin, 


Qidear = o7;'[L sin (y/2)] (17) 


With this, we find an expression for the effectiveness* as 
follows: 


1 
n= Q/Qidea = foe — H)di/sin (y/2) (18) 
0 


n is also a function of the three parameters ¢, y, and 
N.. Results for the fin effectiveness will be presented 
in the following section of the paper. 

Since the use of fins requires that a price be paid in 
terms of weight and material, it is important to know 
the relation of fin thickness to height which leads to the 


* This definition of effectiveness differs from that of references 
2 and 3, where the emissivity € was included in the definition of 


Qideal. 
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maximum heat transfer per unit volume of fin material 
Such a relationship can be derived as follows: Suppose 
that for fixed values of emissivity « and opening angle 


y, the dependence of the fin effectiveness 7 on the con- 
duction parameter NV, has been established. That is, 


n = f(N.) or Q = oT7,'L sin (y/2) f(N.) (19 


For a fixed profile area .1, Lt, we want to find the 


conditions which make Q a maximum. To this end, 


Eq. (19) is rewritten as 


0 = o1,1A,sin (y 2) f(N.) . A,*al,* 
: f ; si kt* 


(20 


Then, to maximize Q, we differentiate 


(5°) “ / 1 df ON, 

= > (21 
Ot /4 f- t dN, Or 
From Eqs. (20) and (21), the value of VV, which yields 


the maximum Q for a given fin volume and fixed « and 
y is found through the condition 


l j l 7 
copt 3 d//dN, 3 On ON, 


To compute the optimum values of N,, it has been 
necessary to numerically differentiate curves of 7 versus 
N.. The results thus obtained will be given later. 

In addition to the overall heat transfer performance 
of the fin as characterized by Q or 7, it is also desirable 
to know which portions of the fin are most effective in 
transferring heat to the environment. To supply this 


information, we may compute the ratio 


eX ex/L 
0 | gdx | (8 — 5C) dé 


S ae 
Ss } gdx | (8 — 3K) dé 
where Q, is the heat loss from the fin surface between 
x = Oand x = x. The contribution of different por- 
tions of the fin to the overall heat loss Q may then be 
shown in the form of curves of Q,'Q as a function of 
x/L. In addition, the local heat transfer g per unit 
area is immediately evident as the slope of these 
curves, since 


gq = dQ,/dx (24 


In order to illuminate the later discussion of results, it 
is also useful to observe that by combining Eqs. (15) 
and (5) an alternate expression for g is found, 


g = eol* — all = eal* — H (15a) 
where a = eforagray body. Eq. (15a) states that the 


local heat flux is the difference between the emitted 
radiation and the absorbed incident radiation. 


Results 


Overall Heat Loss and Fin Effectiveness 
By utilizing the solutions of the governing equations 
(11)—(14) the overall heat loss Q from one face of a fin 
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has been evaluated. The results have been expressed 
as a fin effectiveness n, defined as the ratio of the actual 
heat loss Q to the greatest possible heat loss, which 
occurs when the fin surface is a black radiator having 
infinite thermal conductivity. This information is 
plotted on Figs. 2, 3, and 4 as a function of the con- 
duction parameter \V,. The figures correspond respec- 
tively to emissivity values of « = 1.0, 0.75, and 0.5. 
On each of the figures, there are four curves depicting 
opening angles* of 45°, 60°, 90°, and 120°. 

In Fig. 2 (« = 1) all the curves converge properly to 
the greatest possible heat loss value Q = o7,‘L sin (y/2 
as N, —~ O—1i.e., as thermal conductivity ~ ©. 
However, in Figs. 3 and 4 this ideal heat transfer is not 
reached as V,— 0 because the fin surfaces are nonblack. 

It may be worthwhile to highlight the general trends 
in the results. First, it is seen that the fin effectiveness 
n decreases markedly as the conduction parameter .V, 
increases, and this is true regardless of the emissivity 
or opening angle. Rephrasing this in more concrete 
terms, the curves show that the overall heat transfer 
Q decreases as the thermal conductivity decreases, 
provided that the other parameters are held fixed. 
This reduction in heat transfer is related to the con- 
ductivity-controlled drop in temperature along the fin, 
with a resulting diminution in the radiant emission at 
locations near the fin tip. Another expected trend 
which may be recognized from the figures is that the 
heat loss increases with increasing emissivity; but it is 
important to note that Q is not proportional to the 


emissivity. For example, for V, = 3.0 and y = 90°, 
the heat transfer for « = 1 is 35 percent greater than 
for <¢ = 5. 


By multiplying the ordinates by sin (7/2), it may 
be verified that Q increases with increasing opening 
angle for fixed base temperature and emissivity; but of 
course, the number of fins which may be used on a given 
base surface decreases. For a fixed base surface, the 
heat transferred from m7 fin surfaces arranged with an 
opening angle y/2 is always greater than that from /2 
fin surfaces arranged with opening angle y, provided 
that all fins are of height L. In this connection, it may 
be observed from the figures that for any given V, and 
e, the effectiveness is always greater for smaller opening 
angles. This is explained by the well-known radiation 
cavity effect, whereby the energy emerging from an 
enclosure is augmented by the reflections and _ re- 
reflections which have taken place within the cavity. 
The smaller the angle opening, the greater the number 


of internal reflections. 


Optimizing Conditions 

The gains in heat transfer which the use of fins makes 
possible must be purchased by adding material, and 
consequently, weight. In this connection, it is de- 
sirable to use whatever material is added in the most 
efficient manner. In a previous section of this report, 


* The case of y = 180° corresponds to the situation of non- 
interacting plane fins. Effectiveness values reproduced from 


reference 3 are shown in Fig. 9. 
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it was noted that for a given fin volume and for fixed 
emissivity, opening angle, and thermal conductivity, 
there is an optimum arrangement of fin thickness and 
height which gives a maximum heat loss from the fin. 
The condition from which the optimum configuration 
may be determined is given by Eq. (22). This relation- 
ship has been evaluated by use of the curves of Figs. 
2, 3, and 4, and the optimum values of V, thus com- 
puted have been plotted in Fig. 5 as a function of 
emissivity « and opening angle y. For a given y and 
e, the ordinate provides (V,),,,; and the relationship 
between the fin height Z and half thickness ¢ follows as 


L? = (Ne)opkt/ oT? (25) 


From the figure, the optimum value of N, is seen to 
decrease with increasing opening angle and increasing 
surface emissivity. In concrete terms, this means that 
for a given thickness optimum performance is realized 
for shorter fin heights as the opening angle or the 
emissivity increases. Alternately, if y, «, and /¢ are 
fixed, the optimum height increases with increasing 
thermal conductivity. 

It is interesting tu note that the parameter JV, closely 
resembles the parameter which describes minimum 
weight conditions for a fin which transfers heat from 
its surface to a fluid by convection. With a film heat 
transfer coefficient 4, this parameter, as defined for in- 
stance in reference 11, page 47, reads with the nomen- 
clature used in this report 


u = (L/t)[ht/k]” 
The square of this value, 
u? = L*h/kt 


corresponds closely to V,. The term a7) has the same 
dimension as the heat transfer coefficient h and can be 
interpreted as heat transfer coefficient by radiation of 
a black surface with temperature 7), radiating into a 
surrounding with temperature zero (reference 11, page 
422). The numerical value of «? for a rectangular fin 
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Fic. 7(b) Representative temperature distributions, y = 45°. 


is 2.02 compared to the values 0.9 to 2.3 of the param- 
eter NV, of the radiation fin. The values of both param- 
eters agree, therefore, reasonably closely with each 
other. 

To provide a feeling for the size of the optimum fin, 
let us consider an illustrative example. Suppose that: 
7, = 1000°R (540°F), € = 0.75, thickness = 2¢ = 
0.03 inches, y = 90°, and k = 100 Btu/hr-it-°F. 
Corresponding to the given y and e¢, the value of (V,),,, 
is read as 1.3. Then, from Eq. (25), it follows that the 
optimum height is 3.7 inches. 


Local Heat Transfer 


It is of interest to know the extent to which various 
portions of the fin surface contribute to the overall 
heat transfer. As has already been noted in connection 
with Eq. (15a), the net local heat flux g leaving the 
surface is the difference between the radiant emission 
and the absorbed incidence. For a given fin, it is 
clear from Eq. (15a) that the factors which govern the 
local heat transfer performance are the temperature 7 
and the incident energy /7. Consequently, the vari- 
ation of g along the fin is tied to the variations of 7 
and H. Now, the incident energy is greatest near the 
fin base and decreases monotonically with distance 
toward the fin tip. The temperature is also greatest 
near the fin base and provided that the thermal con- 
ductivity is finite, will decrease monotonically to the 
tip. Consequently, from Eq. (15a), it is seen that g 
may either decrease or increase along the fin (or even 
do both), depending upon whether the 7 or the // 
variation wins out. 

It is possible, however, to outline what might be ex- 
pected for different values of the physical parameters. 
When the thermal conductivity is infinite (V, = 0), 
the temperature will be uniform over the entire fin. 
Consequently, since the incident energy // is largest 
near the base and lowest at the tip, then the heat loss 
is smallest near the base and greatest at the tip. Fur- 
ther, for a given e, the general level of /7 is larger as 
the opening angle grows smaller, and hence, the varia- 
tion of H more greatly influences the local heat transfer. 
Also, for a given opening angle, the magnitude of 
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increases with increasing «. Thus, for the condition 
of infinite thermal conductivity, the greatest variations 
of local heat transfer will occur for small opening 
angles and high emissivities, with g being smallest at 
the base and largest at the tip. For large opening 
angles, HZ is very small and has little effect on the local 
heat transfer. Consequently, g should be almost uni- 
form for large opening angles and infinite thermal 
conductivity. 

As the thermal conductivity is reduced (increasing 
N,), the temperature at the fin tip drops relative to 
that of the fin base, thereby improving the relative 
ability of regions near the base to transfer heat. For 
large opening angles, where /7 plays a small role in 
determining the heat transfer, the effect of temperature 
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Fic. 8. Possible design configurations. 


variation should win out and the heat transfer would 
be expected to have its greatest value at the base and 
least value at the tip. For small opening angles, where 
the level of #7 is large and its variation is of importance, 
gq may either increase or decrease depending upon 
whether the 7 or H variation is stronger. 

These remarks are illustrated numerically in Fig. 6, 
where there is shown a sampling of the local heat- 
transfer results. Space limitations require that the 
presentation be restricted to a few representative situa- 
tions. In the figure is plotted the ratio Q,/Q as a func- 
tion of position x/Z along the fin, where Q, is the heat 
transferred from x = 0 tox = x. The upper part of 
the figure holds for y = 120°, the largest opening angle 
studied, while the lower part is for y = 45°, the smallest 
angle. In each portion of the figure, curves are given 
for «= lt and « = 0.5 fer NV, = 0, 1.25, and’ 40. In 
interpreting the figure, it is to be remembered that 
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along a given curve, the local heat flux is proportional 
to the slope. 

Consider first the results for 120°. 
0), the smallest heat transfer 


For the infinite 
conductivity case (V, = 
does occur in the neighborhood of the base; however, 
the variation of g over the entire fin is very slight. (If 
q = const., Q,/Q = x/L, a straight line with 45° slope.) 
For finite thermal conductivity (V,> 0), it is seen that 
the greatest heat transfer occurs at the base and the 
smallest at the tip. 

We turn next to the results for y = 45°, which ap- 
pear in the lower part of the figure. For the black 
surface (« = 1) with infinite thermal conductivity 
(NV. = QO), the heat transfer in the region adjacent to 
the base is appreciably less than in the tip region. 
The variation of the heat transfer along the fin is not 
quite so drastic for « = 0.5 and infinite conductivity. 
As the thermal conductivity decreases (increasing 
N,), the contribution of the region near the base to the 
total heat transfer becomes relatively larger. How- 
ever, the location at which the greatest heat transfer 


occurs remains at the tip, up to NV, = 1.25, and for 
larger N. values gradually approaches the base—e.g., 
NH, = £8. 


These results for interacting radiating fins are to be 
contrasted with the analytical findings for a convection 
fin, for which the greatest heat transfer is always at the 
base and the least heat transfer is at the tip. How- 
ever, these predictions for the convection fin are based 
on a uniform heat transfer coefficient and uniform bulk 
temperature. If variations of these parameters along 
the fin were included, it is quite possible that results 
would be found which are parallel to these presented 
here for the radiating fin. 


Temperature Distributions 

The variation of temperature along the fin may also 
be of interest. Again, space limitations permit pres- 
entation of only a sampling of the available results. 
Representative temperature distributions have been 
plotted as a function of position x/Z in Figs. 7(a) and 
7(b), respectively for 120° and 45° opening angles. 
In each figure, curves are given for emissivity values 
of 1.0 and 0.5 for NV, = 0, 1.25, and 4.0. 

When the thermal conductivity is infinite (V, = 0) 
the temperature is uniform along the fin regardless of 
finite thermal 
decreases 


emissivity. For 
temperature 


opening angle and 
conductivity (V, > 0), the 
monotonically along the fin; and as the conductivity 1s 
reduced, the temperature variation is accentuated. 
For a fixed value of NV, the heat loss from the fin in- 
creases with increasing opening angle and increasing 
emissivity. Consequently, the temperature drop be- 
tween the base and the tip is magnified as these param- 
eters increase. Numerical verification of these remarks 
may be found in Figs. 7(a) and 7(b). 

The drop in temperature from the fin base to the fin 
tip is magnified in its effect on the radiant emissive 
power by the fourth-power relationship. For a tem- 
perature 7/7), = 0.65 at the fin tip, the emissive power 


OF 
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at that location is only 18 percent of that at the base. 


Applications 


Fins on a Plane Wall 

The analysis has been carried out with primary em- 
phasis on the situation of a symmetrical heat transfer 
from both faces of a fin. For applications in which one 
face of the fin experiences a much greater heat loss than 
the other, the same analysis still applies, with the excep- 
tion that ¢ is now interpreted as the total fin thickness. 

Now, consider an ensemble of longitudinal fins erected 
on a plane wall, as shown in Fig. l(c). In the absence 
of strong external radiation it is to be expected that the 
heat transferred from the fin surface which faces the 
outside environment will be much greater than that 
from the surface which sees the plane wall. So, the 
results presented in this paper apply approximately. 
In the introduction, certain qualitative remarks were 
made about the conditions under which the use of fins 
would either increase or decrease the heat transfer rela- 
Now, we are in a position 
The heat transferred 


tive to the unfinned base. 
to give quantitative information. 
from a section of unfinned base surface of dimensions 
L sin (y/2) X lis 
eo1,'L sin (y/2) 

This is to be compared with the heat Q transferred from 
a fin of length Z arranged at an angle y/2 with the 
vertical. Thus, fins decrease the heat 
transfer depending upon whether 


increase or 


Q/eaTy'L sin (y/2) 

But, the effectiveness 7 has 
previously been defined as » = Q/o7,'L sin (y/2). 
And so, the criterion for comparing the heat trans- 
ferred from the finned surface with that for the un- 


is greater or less than one. 


finned base is 


7 > & 


If 7 > &, then the finned surface transfers more heat 
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(reference 3). 
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than the unfinned base; while if 7 < «, finning actually 
reduces the heat transfer. The break-even point is 
Nn = &. 

With this, we can return to Figs. 2, 3, and 4 and find 
a quantitative measure of the conditions under which 
finning has utility in increasing the heat transferred. 
For simplicity, suppose that the fin and base surfaces 
have the same emissivity. First, for «e = 1 (Fig. 2), it 
is seen that at best, 7 < «. Consequently, finning 
will always diminish the heat transfer except for the 
condition of infinite thermal conductivity (V, = 0). 
Next, for « = 0.75 (Fig. 3), the break-even point ranges 
from V, = 0.09 to 0.29 as the opening angle ranges from 
120° to 45°. For N, values less than the break-even 
point, finning is advantageous; while for N, greater 
than the break-even value, finning decreases the heat 
transfer. So, for example, for the opening angle of 
15°, fins increase the heat transfer only when 0 < N, < 
0.29. Finally, turning to Fig. 4 for e« = 0.5, it is seen 
that the break-even point ranges from N, = 0.12 to 
1.06 as N. ranges from 120° to 45°. The conclusion 
to be drawn is that for the condition where the fins 
have the same emissivity as the base surface, there is 
little possibility of surpassing the performance of the 
unfinned base except where « is relatively low. How- 
ever, if the base surface has low emissivity while the 
fins are made of highly emissive materials (« > ¢,), 
then the use of fins may significantly increase the heat 
transfer. 

Design Considerations 

The analysis as presented in this paper has assumed 
that the finned surface can radiate freely into a sur- 
rounding of low temperature. In reality other objects 
may be located nearby which have a temperature suffi- 
ciently high that the radiation received by the finned 
surface from these objects cannot be neglected. A 
similar effect will result when a heat exchanger consists 
of several elements of finned surfaces. The design of 
such a heat exchanger has to take into account the 
mutual shading of the finned elements and keep it small 
by proper arrangement. One configuration which 
suggests itself for a heat exchanger consists of finned 
tubes with the tubes arranged in one plane and parallel 
to each other and interconnected by fins which are 
either arranged in the same plane or located in different 
planes. Two configurations which can be analyzed in 
a simple way by application of the results of this paper 
are shown in Figs. S(a) and S(b). 

For these two arrangements, the projected lengths 
P and P, (in the plane of the tubes), and the fin weights 
IV and W,, respectively, may be compared, assuming 
that the two heat exchangers transfer the same amount 
of heat Q with the same base temperature 7), emissivity 
e, and conductivity k. The fin height Z may also have 
a fixed value. The thickness of the crimped fin may 
be ¢ and the thickness of the straight fin 2¢. The 
amount of heat transferred by a projected length P and 
unit width of the heat exchanger in Fig. 8(b) is 


Q = nol;'P 
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Correspondingly, for the exchanger with straight fins 
Q = npol,'P> 


The effectiveness n, of the exchanger in Fig. S(a) has 
been calculated in reference 3 and is presented in Fig. 
9. From the analysis it follows that one curve is 
sufficient to represent the effectiveness of this geometry 
for any « when the parameters as indicated on the 
ordinate and abscissa of Fig. 9 are used. From the 
above two equations, the ratio of projected lengths is 


obtained as 
P P, = %/9 (26) 


Introduction of numerical values for n from Figs. 2—4 
and of n, for the same value of V, from Fig. 9 demon- 
strates that the projected length of an exchanger can 
be decreased by crimping of the fins. This may be 
important for structural considerations or in the effort 
to shield the exchanger from solar irradiation. The 
fin weight of the exchanger with plane fins is W, 
2tP,d (d = density of the fin material). The heat 
flux Q can, therefore, be expressed by 


Q = npol,*(W, 2dt) 
The weight of the exchanger with crimped fins is W = 
2¢Pd/sin (y 2) and the heat flux is 


W sin (y/2) 
) = Co . 
‘ ? 2td 


From these two equations, one obtains for the ratio of 
fin weights 

W/W, = np/n sin (y/2) (27) 
By insertion of numerical values into the right-hand 
side of the equation, one finds that the fin weight of the 
exchanger with crimped fins is always larger than the 
fin weight of the exchanger with plane fins. This 
indicates that one has to pay a weight penalty for the 
fact that the projected area of the exchanger with 
crimped fins is smaller. The decision as to which of 
the two features is the more important one will depend 
on the specific situation. The fact that Fig. 8(b) 
provides some protection of the tubes against the impact 
of meteorites may count in favor of the exchanger with 


crimped fins. 
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Bending of Circular Plates Due to 
Asymmetric Temperature Distribution 


MARVIN FORRAY* ann MALCOLM NEWMAN** 


Republic Aviation Corporation 


Summary 


This paper considers the problem of a circular plate sub- 
iected to a thermal gradient which varies linearly through the 
thickness and is asymmetrical over the planform. The general 
Eqs. (1)-(9) presented here are applicable both to the solid circu- 
lar plate and to plates having concentric holes. A detailed solu- 
tion and nondimensional design curves are given for the clamped 


solid circular plate. 


Symbols 
radius of solid plate 
D plate flexural rigidity, Eh*/12(1 — v? 
E = Young’s modulus 
h thickness of plate 
km = nonnegative integers 
V,,Mo bending moments 
Me = twisting moments 
q = normal load intensity 
Q,,Qg = transverse shearing forces 
r, 0, 2 = cylindrical coordinates 
t = time 
Tp = temperature difference between upper and lower 
faces of plate 
V; reduced shear 
ti = deflection 
a = linear coefficient of thermal expansion 
v = Poisson’s ratio 
Cj stress components 
Introduction 


Pigeon PLATE-LIKE BULKHEADS in air and space 
vehicles are subject to temperature distributions 
which are rotationally unsymmetric and vary through 
the thickness. The temperature distribution in such a 
plate is expressible in the form of a power series in the 


thickness coordinate—1.e., 
i = > T(r, 0,t)2" 
k=0 


For thin plates the temperature may be approximated 


by 


T(r, 0,t) + 27\(r, 0,t) 
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NS) 


where the effects of the two terms may be superposed 
if the deflections are small and the material is linearly 
elastic.’ The first term, 7)(7, 0,/) gives rise to a plane 
stress or slab problem. Solutions for the quasi-static 
case of an elastic plate with uniform thickness subject 
to an asymmetric temperature variation 7)(7, 96,0) 
over the planform, were considered by Horvay? for 
solid circular plates and Forray* for rings. The second 
term 7}(r, 0,t)z gives rise to a bending problem. It is 
the purpose of this paper to consider the deflections 
and stresses in circular plates due to such a linear ther- 
mal gradient through the thickness when the tempera- 
ture distribution is arbitrary over the planform —i.e., 


T(r, 6,to)2 = (—2/h)Tpl(r, 0 


where / is the thickness of the plate and 7» is the 
temperature difference between the upper and lower 
faces of the plate (Fig. 1). The assumptions of classical 
thin-plate theory will be employed. ! 

It will be shown that the same governing differential 
equation applies to the plate bending problem as to 
the plane stress problem. A detailed solution and 
design curves will be given for the clamped circular 
plate. The authors are at present developing the solu- 
tions corresponding to other boundary conditions in 
both solid and hollow plates. 


Plate Equilibrium Relations 


In cylindrical coordinates, the stress equilibrium 


equations are‘ 


Oo;, l Odor Oo;- Orr — Coe 
+ + 0 
or r O86 Oz r 
Ocre ‘ 1 Oag O09: + 2or9 " 68 
Or r O80 Oz r 


Oo;; 4 1 Oop; 4 Oo.- G; 
Or r O86 Oz r 


The plate (Fig. 1) is assumed to be subjected to tem- 
perature and distributed normal forces over its faces 
The first two equations of (1) are multiplied by z, then 
integrated through the thickness, while the third equa- 
tion is integrated through the thickness, yielding the 
plate equilibrium equations 
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ON, 1 OM, + M, — Ms, - 
Or r oO ¢, r - 
ON, 10M, 2 
— + —-Q— My = 0 (2 
or r O86 ‘ r , | \<) 


O00, 100, . @, 
‘ = (0) 
or J r OO " r 7? 


where g = 6:22\2=n/2 — Gz2| 2 
shears per unit of length are defined by 


h/2 2H, 2 
M, = zo,, dz M, = | O99 dz 


h/2 « 


h/2 h/2 
0, = f Or, dz Op = | C2 dz 
—h/2 —h/2 | 


) 


From the assumptions of classical plate theory’ 

(1) The plate is assumed to be homogeneous and 
isotropic with temperature-independent material prop- 
erties: 

(2) Plane sections which before bending are normal 
to the median plane of the plate remain plane and nor- 
mal to the median plane after bending; 

(3) Hooke’s law holds, so that the bending and twist- 
ing moments per unit length can be expressed in terms 
of the normal deflection w and the temperature differ- 
ence 7'p(r,@) between the upper and lower faces. 

Substitution of the known expressions for M7/,, Mg, 
M,, into the first two equilibrium Eqs. (2) yields Q, and 
Q, respectively. Thus, 


M, = -—D E +? (‘% + =e - 
r a 
aTp(1 £7) 
h 
M, «= <9 E oe +i ~ 
i y= 


aTp(l + | 
—— - (4) 


h 


My = D(1 — v)(0/0r)(we/r) 
l Tp 
i oaie | v2 +e | 
or h 
1O0 (1 + v)aTp 
» = —D- =| Vw —- ———— 
Co r al h | 


Kirchhoff’s reduced shear V, = Q, — (1/r)(0.M,9/08) 
is given by 


U'ee 


(1 — »)Dwep 


: (4a) 
yt 


(i + pel ™ 
h 


y»/2,and the moments and 
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Fic. 1. Plate geometry. 


When the expressions for Q, and Q, from (4) are sub 
stituted into the third of Eqs. (2) there results (with 
q = VU) 


l+p = 


Viw = V2(aTlp) (5) 
h 


where V? is the two-dimensional Laplacian operator 
and A* = V°V*. This equation is of the same form as 
the equation for the Airy stress function in the plane 


stress problem.°® 


Solution of the Deflection Equilibrium Equation 
It is assumed that (1 + v)aZ7‘p/h is expressible in the 


form 


l+p 
aTp = LL DX Aum r* cos mé + 


h m=0 k=0 


> > 2; Bynr® sin m6 (6) 


m=i k=0 
For the general term 


(1 + vial p . 
— = A,» * cos m6 + Bin, r* sin mé 


h 
the deflection w may be expressed as the sum of the 
complete solution of V‘w, = 0 and a particular solution 
of V?w, = (1 + v)alp/h.*® The expression for w, is 


We = Ayo + bor? + cor? In r + dy In r + 
(ayr + br? + cr! + dr in r) cos 6 + 
(ay’r + by’r? + cy’r-! + dy’r In r) sin @ + 


oo 


2. [(a,r” -+ br" t+? + c,r-" + d,r—"*?) cos nO + 


n=2 
(a,/r" +: b,'r" +2 + c,'r-" + d,’r-"*?) sin nO] (7) 


If the boundary conditions are axisymmetric, only 
the antisymmetric terms in w, will occur for temperature 
differences in sin m6, while the symmetric terms of w, 
will be associated with temperature differences in 
cos mé. Consequently, for 

_ ( Rn 
(] a val D | Ate ; Coe me 
h Bim, r* sin mé 








sub 
vith 


itor 
l as 
ane 


on 


the 


(6) 


he 
on 





BENDING OF 


° ° ° P 25 
the particular solution is given by 


2m (r) cos mé 


Wy ; (S) 
F h,,(r) sin mé 
where 
. . | 
7 2 
Bm _ y—m | y2m | aw gk t+l—™ dr dr 
hy, Bum 
The integration yields 
Apm rt? cos m0 
(k + 2)? — m* 
W, = (9) 
B,.,, r**? sin mé 
(k + 2)? — m? 


when k + 2 — m # 0), and 


07 r T T T m=O [> 


p*+2 


Axo 











.006 m=2 


004 4 


002 





-.002 


Ww 
Axe b**® cos 26 





-.004 
yp kel 


-.006 








-.008 
O z 4 6 8 10 


r/b- 


CIRCULAR PLATES why 


In + l 

I kes a - |r” cos mé 
2m (2m)° 
In 7 l ; 

Bim _ | r™ sin mé 
2m (2m)- 


when k +2—m=0. The deflection w = w, + w, 


W, = 
p (Qa) 


Clamped Solid Circular Plate 


For a solid plate, the constants c,, c." dy a 
0, Eq. (7), from the requirement that the deflection and 
moments remain finite atr = 0. The constants a, and 
b, are determined from the appropriate boundary condi 
tions. For the clamped plate (w|,.» = dw/dr|,» = 0 
the expressions for the deflection, moments, and shears 
in nondimensional form are as follows 








014 me=l 
012 5 
—4 
010 3 
@ 2 
3 008 
3/3 
oe 
2 
= .006 
< 
004 
002 
K=| 
fe) 
fe) 2 4 6 € 1.0 
/b 





-.002 


w 
Axsb"*®cos 36 


-.004 





-.006 











-.008 
.¢) 2 4 & 8 1.0 


r/b 


Fic. 2. Nondimensional deflections as a function of distance from the center of plate 
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BENDING OF CIRCULAR PLATES ivi 


Case .A: = 1p, r°cosmé, k+2—m# (0 


ti! ; | (‘ — m\ [r\" k+2—m r\"+2 r\k+2 
Limb t? (+ 2)? — m? 2 b 2 b aa aained 
Af, ; | k—m ; de i ke+2—-—m 
| DB . aa an —m(m — 1) 5 (1 — pv) b + (m + 1 9 x 
] ») »))\t r E , {r 
(8 -- 5) WM Zi i +(1—v(kR+2—m ¥ cos mé (10) 
M, | k—m r\""> + (m +1) (kR +2 —m 
; mim — 1) (l — » 
PS a (k + 2)? — m- 2 b 2? 
) » J +)! r\" r\ 
»(m — 2) + vim + 2), + (1 — v)(kR + 2)(kR+ 1) {- cos mé 
b \0 
M,e —({l—~y | ' (F —m (: m— : (* +2—m (; " 
= mm — ) — m(m ) + 
LimDbY — (k + 2)? — m? 2 b sl 2 b 
r\' 
m(k + 1) ( ) | sin mé 
b 


0, 2m(m + 1)\(k + 2 —m) [r\"" 
1 cos me 
AreD (k + 2)? — m* b 
Oo —2m(m + 1)(k +2 —m) [r\""' . , 
in mé 
A pmb (k + 2 -_— m* b ° 
The reduced shear at the boundary is 
ai 2m(m + | , 
= cos mé (10a) 
ADF ke+2+m™ 
: (1 + val p 
Case B: > _l,, r* cosmé, kR+2—m () 
h 
WW l j r | r m r m+2 
- = 1+ 21n 4 — cos mé 
Apmb§ +2, Am LA b/§ \b b 
M, | — (: ae 1) Se +2 ») (’ 
= —(m — 2)(m + 1)(1 — »v) + (m+ 1) 4(m 2) — vim — 2 -- 
NicmDb® Am b '\o 


r gl 
(Q2m)(m — 1)(1 — v) {In F F cos mé 
) ) 


) +m +1) {=m — 2) + v(m + 2)} (Z) + 
) 


r 7r\"—2 
(2m)(m — 1)(1 — pv) (in )( cosmé (11) 
b, b 
Myo (1 — v) (m+ 1) "é We j 2m — 1). ry fr\™? 1. 
= — 41+ in -7 sin mé 
A pmDb | b (m + 1) bf \b 


QO Sa 
= (m + 1) cos mé 
ArmDb*—" b 


Og r m 1 ; 
| Dbe= = —(m + 1) b sin mé 
<1 ptm ) 


r 


b 


Mo l 


(m= + 3m — 2)(1 — »v) 
A,,,.D6' 4m | . ( 





The reduced shear at the boundary is given by 


Vs r=b 
= (m + 1) cos mé (11a) 
A pmDB*—! 
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Fic. 4. 


Note: 
by sin m6, then for the nondimensional w, M,, M,, 
Q,, and V,, replace cos m@ by sin m6, while for Qs, 
M,9 replace sin m@ by —cos mé. 

The results of Eqs. (10) and (11) are given in Figs. 
24fork = 1,2 ..... 5 with m = 1,2:3 and for & = 0,1,2 

. ) with m = 0 (axisymmetric case of reference 7). 


If, in the expression for 7'p, cos m@ is replaced 
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Longitudinal Dynamics of a Lifting Vehicle 


in Orbital Flight' 


BERNARD ETKIN* 
Institute of Aerophysics, Unwersity of Toronto 


Summary 


\ report is presented of an analytical study of the longitudinal 
perturbations to the steady orbital flight of a lifting self-propelled 
vehicle in a circular orbit. The equations are so formulated as to 
lead to the correct limits at the two extremes 
Additions to the conven- 


satellites in space 
and aircraft in the lower atmosphere 
tional equations of airplane dynamics are found to be necessary 
to account for the vertical gradients of gravity and atmospheric 
density and for the curvature of the undisturbed flight path. 
Calculations were performed with both linear and nonlinear equa 
tions for a particular set of aerodynamic and inertial parameters, 
and altitudes from 100,000 to 700,000 ft. The linear solutions 
contain two oscillatory components which can be identified with 
the classical phugoid and short-period modes, and a new spiral 
mode which may be unstable or stable, depending on the thrust 
law. With increasing altitude, the period of the phugoid ap- 
proaches asymptotically to the orbital period, and that of the 
“short’’-period mode tends to infinity. In the altitude range of 
nearly equal periods there is beating of the two oscillations, pro 
ducing dynamic instability in pitch. The solutions of the non- 
linear equations show that the linear theory is valid for a practical 


range of the disturbance amplitudes 


Symbols 
A, B, C = moments of inertia 
c = mean chord 
Cp = drag coefficient 
CL = lift coefficient 
ca = pitching moment coefficient 
Cr = thrust coefficient 
Cw = mgo/(1/2)pouo?2.S 
C,, C. = coefficients of aerodynamic force in directions of 
body-fixed axes 
Cn = OC,,/00 
Cue = OC,,/d4 
Cones = 0C,,/d0a 
Cné = 0Cm/0(ae/2u0) 
Con = 0C,/da 
Coe = 0C,/da 
Cra = oc, /Oa 
D = drag, or d/di 
gy = acceleration due to gravity at reference altitude 
g’ = perturbation to acceleration of gravity 
ig = yu(k,/1)? 


Presented at the Aerodynamics Session, IAS National Summer 
Meeting, Los Angeles, June 28—July 1, 1960. Revised and re 
ceived October 10, 1960. 
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Research Board of Canada is also gratefully acknowledged. 

This paper consists of a condensation and extension of ref- 
erence 8. The analog computer results are due to R. S. Rangi.'° 

* Professor of Aeronautical Engineering 


k, = radius of gyration in pitch 

L = lift 

l = characteristic length of vehicle, €/2 

m = mass of vehicle 

M = pitching moment 

Q = angular velocity in pitch relative to the earth 

q = perturbation of Q 

q = gl/2u 

R = position vector of vehicle 

ro = radius of reference circular orbit 

’ = perturbation of R 

S pags 

r = 7/A 

Ss = speed ratio, uo/m, 

= = reference area 

2 = thrust 

t = time, sec. 

s* time unit, //™ 

i = ¢/t* 

U,W = components of v’c in the directions of the body 
fixed axes 

u, Ww = perturbations of (U, W) 

Ug = reference flight speed 

Uc = circular speed [Eq. (2.3)] 

a = u/uy 

vc = velocity relative to earth-fixed axes 

oe 4 = aerodynamic forces in directions of body-fixed axes 

“, 2 = body-fixed axes 

a = W/uUo 

B = (1/p)(Op/or) 

8 = angle of pitch (see Fig. 5) 

0 = perturbation of 0 

w = angular velocity of the earth 

€ = m/pSro 

p = air density 

Po = air density at the reference altitude 

mn = m/pSl 

A = wavelength of flight path in an oscillatory normal 
mode 

d = root of characteristic equation 


(1) Introduction 


A NUMBER OF ANALYTICAL STUDIES have been pub- 
lished on the dynamics of flight at high speed and 
high altitude—e.g., references 14. These deal with 
transient, unpowered flights of the re-entry or boost- 
glide type, and have done much to provide an under- 
standing of flight in this régime. 

However, it was felt that worth while 
understanding could be derived from another ap- 
proach—i.e., from an extension of the classical analysis 
of small disturbances of steady powered flight up to the 
speed and altitude at which satellites operate. The re- 
sults of such an analysis would be to some extent rele- 


additional 
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Fic. 1. Thrust required for steady flight in circular orbit. 


vant to the shallow entry of a lifting vehicle, or to a 
boost-glide vehicle in the glide phase, in addition to 
being directly applicable to the sustained flight of a 
self-propelled sub-satellite describing a steady circular 
orbit. 

Such an investigation is reported here. It consists of 
(1) the calculation of the properties of the normal modes 
of a particular configuration from 100,000 to 700,000 ft 
altitude at a small constant lift coefficient, and (2) the 
calculation of transients with second-order terms in- 
cluded. The effects of the thrust law upon the modes 
are examined, by considering the cases 

(1) constant-thrust rocket engine 

(2) air-breathing engine (thrust « density) 

The thrust which the engine must supply becomes 
very small at high altitude, as shown in Fig. 1. At 
350,000 ft, the thrust weight ratio is already down to 
ig. 

In order to obtain a set of equations which give cor- 
rect results in the upper atmosphere, and which, as well, 
lead to the correct limits at the two extremes-—-satellite 
in space and aircraft in the lower atmosphere-—it was 
found necessary to make some additions to the conven- 
tional equations of airplane dynamics, as given for 
example in reference 6. These additions account for the 
gradients of density and gravity with height, and for the 
curvature of the steady reference flight path. Neu- 
mark* had shown previously the necessity for including 
the density gradient. 


(2) The Equations of Motion 


The equations of motion used here apply strictly to 
a nonrotating spherical earth, with an atmosphere at 


SCIENCES—OCTOBER 1961 


rest. However, the effects of the rotation of the earth 
are not very important. Fig. 2 shows that the Coriolis 
component of the acceleration remains less than 10 per- 
cent of g up to satellite speed; and the centripetal ac- 
celeration due to earth rotation can be absorbed into the 
local value of g without introducing appreciable error 
The effects of relative motion of the atmosphere could 
be treated as gustlike disturbances to an otherwise 
steady flight. Such data as exist indicate that mean 
winds to not exceed a few hundred feet per second at al- 
titudes up to about 100,000 ft——1.e., they are small com- 
pared with the vehicle velocities envisaged 


The Reference Steady State 
For steady motion on a circular path (Fig. 3) the 
equations of motion become 


Muy”, et 


M80 _ Lo . (2 1) 
M, = Of 


Ty = Do 
where the subscript zero identifies the reference condi- 
tion. 

From the first of Eqs. (2.1), 


Lo = mgo(1 — mo? rogo) (2.2) 

When L, = VU. we have the condition of circular orbital 
flight—.e., 

TP V rg (2.3) 


We define the ‘speed ratio” as 


S = U/l, 
so that Lo = mg(1 — s?) (2.4) 
whence C,, = 2e1 — s*)/s? (2.5) 
where e = m/pSro = pl/n; 


and yw is the customary relative-density parameter 
m/pSi. Eq. (2.5) is plotted on Fig. 4 with ¢ as the 
ordinate and the speed ratio s as the abscissa. For all 
values of C;,, ¢€— © ass— 1.0. For any given wing 
loading, and atmospheric structure, the scale of «€ be- 








Fic. 2. Acceleration components, percent of ¢ 
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Fic. 3. Reference steady state 


comes a scale of height. This is shown at the right in 
Fig. 4 for the ARDC model atmosphere’ and W/S = 50 
psf. It is seen that for all practical lift coefficients the 
speed is very nearly equal to uv, for heights above 300,- 
000 ft. Between 100,000 and 300,000 ft is the zone of 
most rapid change in s 


The Small-Disturbance Equations 

The equations of motion are written in scalar com- 
ponents referred to body axes fixed in the vehicle (Fig. 
5). Considering only motion in the vertical plane, we 
obtain the following equations: 


Xo + AX — mg sin 9 = m(U + QW) 
Zo + AZ + mg cos 8 = m(IW — QU) (2.6) 
My, + AM = BOQ 


where the subscript zero indicates the reference condi- 
tion, and the prefix A indicates the perturbation. !n 
the reference state we have also 


6@=6=0, W=~y 0, R=* 
Q=q, U=u, g= x 
With these values inserted, Eqs. (2.6) yield 
X, = 0 
Zo + mgo = —MUQo (2.7) 
M, = 0 


With these results, Xo, Z,, 9 can be eliminated from 


Eqs. (2.6) to give 


AZ + m(go + g’) cos 0 — mg 
= m[w — (go + g)(uo + u)| + is ‘* 
AM = Bg 


AX — m(go + g’) sin @ mu + (go + q)w] | 
he 


8) 


We now linearize Eqs. (2.8). In so doing we note that 
go is of order 10~* rad. /sec. This is a very small angular 
velocity, and hence it might be expected that its product 
with the small-disturbance quantities u and w would be 
negligibly small. In fact, this has been found to be not 
so in the high altitude limit, and the terms referred to 
must be retained. This implies that the terms gu and 
gw might need to be kept for values of g of the same 


order as go—i.e., that the linearized solution is valid 
only for g < qo, at least above a certain minimum height 
This implication has been investigated by studying 
transient solutions on an analog computer [Section 
Shi. 
On using the inverse square law for g, we get g’ 

—2(go/ro)r, and after dropping all second-order per 
turbation terms, we get from Eq. (2.8) 


AX — mgd = m(u + gow) 
AZ — 2mg(r/ro) = m(w — qou — uUoq) (2.9 
AM = Bg 
The system of equations is completed by the addition of 
the following kinematic relations (see Fig. 5): 
R = Usin® — Wcos® 
which for small disturbances becomes 
r= uf — w 
U cos 80 + W sin 9} 2.10 


and )=@e- 
‘ R 


which when linearized reduces to 
q = 6 — (u/ro) — go(r/ro) (2.11 


The expressions for the force and moment perturba- 
tions, Eqs. (2.12), contain the usual aerodynamic 
derivatives and, in addition, the derivatives Y,, Z,, M, 
which are associated with the vertical gradient in den- 
sity, and the moment derivative WV, which results from 
the gravity gradient [Eqs. (2.12) |: 
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| 
Vf : ash g 
t 
b 
Fic 5. Body axes and notation. 
| M 
A 
. - , r sl 
AX X,u + X,w + X,r } 
AZ = 24+ 240 + 27 (2:72) 
~ ee ‘ Fic. 6. The dumbell effect: the gravity gradient 0g/0r produces 
AM = Maw + M.w + Mg + Mx + Me) none. ee 
When Egs. (2.12) are introduced into (2.9), the com- b 
plete set of linear equations, in matrix form, is 
d y y é a ae T 
m 7 = Ae (mgo — Aw) 0 Mg —X, Uu “a 
. d ' 2mgo vi 
—(Z, + mq) m — Zu —mMuy 0) —Z, + w 
dt ro 
0 ( nt oe (8 «8 ul ul ; 
— {| M; M = i. =—% —M,;, ¢ (0) 24a) 
* dt dt : : / (2.18) 
l d aT 
Q l — t a] 
To dt a 
d 
0 l 0 — Uo r 
t dt | Ll J 
Values of the r Derivatives 
M,. The aerodynamic pitching moment is given by 
M = C01/2)."5€ - 
and therefore (OM/Or)) = C,,(1/2)v,2SE Op/Or va 
= C,,,(1/2)uo2SE Op/Or 
Since C,, = 0 in the reference condition, then /, = 0. 
Z,. Similarly, 
Th 
Z, = C,,(1/2) potto2S(1/ po)(Op/Or) 
The bracketed term is denoted 8, and C,, = —C,,, so that th 
_ wii shi 
Z, = —BC,,(1/2) pouo?S (2.14) ie 
X,. The X force is the difference between the thrust and the drag, ert 
, am , on qu 
As fT —D and os T= D th 
Calculating as for Z, above, the value of D, is tiv 
tio 
D, c= BCp,( 1/ 2) potto2.S 








luces 


14) 
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For the thrust we consider two possibilities: 


|) rocket engine of constant thrust; 7’, = 0 and 


xX, = —D, = —BCp,(1 2) potto”.S 


(2) air-breathing engine with 7° « p, so that 


- = BCr,(1 2) potty?.S 


Since Cp, = Cp,, then for this case 
X, = 0 
\J,. Because of the inverse square law, mass elements which are farther from the earth are subject to a smaller 


When the axis of the dumbell is horizontal 


gravitational attraction. This is illustrated for a dumbell in Fig. 6. 
For an arbitrary 


the equilibrium is unstable (./, > 0), whereas the vertical orientation is a stable one (VW, < 0). 
body, the expression for ./ is® 
M = —(3/2)(g-/R.)(A — C) sin 2x 


where g,, RX, are respectively the values at the mass center of the gravitational force per unit mass, and the radius; 
and x is the angle of elevation of C,,. It is as- 


A, C are the moments of inertia about the principal axes C,,, Cz,; 
sumed here that the x, z axes are principal axes, and hence that x = 9. For small angles, then, 
M = —3(g0/r)(A — C)@ and Me = —3(g0/r70)(A — C) (2.16) 
If the body in question were a uniform slender rod (C < 1) oriented vertically, the period of small oscillations would 
be 
T = 24V A/(—M,) = (24/V3)V ro/go = (1/03) T orbit 
Thus, although it is very small, the gravity-gradient effect is nevertheless capable of producing technically important 


oscillations. At great enough heights it becomes more important than the aerodynamic pitching moment. 


The Nondimensional Equations 
When the expressions for the derivatives are substituted into Eqs. (2.13), and the latter are then made nondimen- 


sional by the definitions in the list of symbols, they become 


, ' . Cp,8r i | 
(uD — Cx) ~ (26 + Cua 0 Cw vos |i? 
2(Cz, + €) (2uD — X,,) —2u 0 (2Cw + C.,Bro| | a 
0 —(CreD + C,,) (tpD — C,,,) —Cm, 0 i. 0 (2.17) 
/ ! D 
ro ro 
f 0 l —D - , la 
ro Yo L 


The alternative values for the last element of the first row correspond to the two alternatives for the thrust, the zero 


value being for the air-breathing engine. 
Ea. (2.17) is the system that was used for the calculations of the small-disturbance modes. 


(3) Numerical Data Used 


The Atmosphere ’ 
pa , z : not reproduced here. However, the values of 8, which 
The ARDC Model atmosphere’ was assumed for all é : ive * 

‘ ‘ ta are not given there, were calculated from the tables of p 

the calculations. Available test data from probe firings ie a aa ' é' 

sg ' = : and are plotted in Fig. 7. The kinks in the curve occur 
show that it is quite good up to about 300,000 ft, but Pie 
; : : 6 oa ‘ at altitudes where there is a kink in the assumed tem- 
deductions from satellite flights indicate it may be in fi {Hai 
; ? : perature profile, and are presumably not very realistic. 
error at much higher altitudes. In any event, only the 
quantitative details of the results given here depend on The Vehicle 
the particular atmospheric model used. The qualita- on . : ; — 
. ; rhe relevant geometrical and inertial parameters as- 
tive results would not be affected by a different assump- 
5 sumed are: 
tion. 
7 L = (1/2)é = 25 ft 


Details of the atmosphere given in reference 7 are 
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Fic. 7. Density-gradient parameter. ARDC model atmosphere. 


B = (1/p) (Op/Oh). 


es = 25 ft. 
W/S = 30 psf (at sea level) 


From these, we get 
uw = 0.0373/p 


The stability derivatives used are listed below. They 
are representative of values obtained by the simple 
Newtonian impact theory for a slender body (cone or 
wedge of about 3° semiangle) at moderate angle of at- 
tack. The exact values used are not in themselves im- 
portant. The main results of the paper would be in- 
sensitive to moderate changes in the numbers. 


Cire = —Cz, = 0.329 a = 0.152 rad. 

Cp = 0.0019 + 0.493a? C,, = 0.05 

Cu, = —0.028 Cp, = 0.0133 

Cm, = —0.0548 Cr, = —2Cp, = —0.0266 
Cm, = 0 Cr, = Cy, — (0CD/da) 
Cmg = 1.41Cywl/r0 = —0.100 


The values of the main parameters that change much 
with altitude are given in Table 1 for heights from 100,- 
000 to 700,000 ft. 


(4) Normal Modes 


The eigenvalues and normal modes (eigenvectors) of 
the system of Eqs. (2.17) were found with the use of an 
IBM 704 computer.* As noted in Section (3), the 
value of C, was held constant at 0.05. Four different 
cases were calculated: 

(Aj Complete equations, rocket engine (/’ = const.) 


(B) Complete equations, air-breathing engine 
(T x p) 

(C) Approximate equations with 8 = 0 (same for 
both engines) 

(D) Approximate equations with g = 0 in Eq. 


(2.13), rocket engine. 
In each case, calculations were carried out for a range 
of heights from 100,000 to 700,000 ft. 
Case (A) provides a set of basic solutions which are 


*The program used was devised by G. Jacquemin and J 
Galipeau, of AVRO Aircraft, to whom the author is indebted for 
valuable assistance. 
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used as a reference for comparison; case (B) shows the 


effect of the thrust law; case (C) shows the importance 
of the density gradient term; and case (D) shows the 
effect of the small gq, on the solutions. 

The characteristic equation of the system is of the 
fifth degree, leading usually to three normal modes, two 
oscillatory and one nonoscillatory. The two oscillatory 
modes can be identified with the classical phugoid and 
short-period modes. The latter is a misnomer at high 
altitude, since the period tends to infinity, exceeding 
that of the phugoid. A better name would be the 
“attitude” or “pitching’’ mode. The nonoscillatory 
mode describes a spiral flight path, which may be stable 
or unstable, depending on thrust law and speed. Up to 
about 400,000 ft altitude it and the phugoid are ‘‘tra- 
jectory”’ or “arrow” modes, in that the angle-of-attack 
perturbation is negligible —i.e., the vehicle flies with its 
axis aligned along the flight path. The main charac- 
teristics of the three modes are given in Figs. 8-12. 
Periods of Oscillatory Modes 

Fig. 8 shows the periods of the two oscillations, and 
how they increase with altitude. The results for the 
short-period mode for cases (A), (B), and (C) were the 
same. Furthermore, the same curve was obtained from 


the equation 


r= 2t* ( = yo (4.1) 
tin, + Cu) 


This is the approximation based on a single degree of 
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Fic. 8. Period of longitudinal oscillations. Cro = 0.05 
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freedom, as though the vehicle were constrained only to 
rotate in pitch about an axis through the mass center. 
The period is seen to increase rapidly with height, cross- 
ing over the phugoid period (or 7orit) at 490,000 ft and 
tending asymptotically to infinity at 505,000 ft. Above 
this altitude, the unstable gravity-gradient effect over- 
powers the aerodynamic pitching moment, and the 
pitching mode becomes divergent. Altitude stabiliza- 
tion in this condition would have to be provided by 
other than aerodynamic means. 

The limiting behavior of this mode at high altitude 
is seen to be governed by the gravity effect. In the ex- 
ample here, C,,¢ is positive (unstable) and the period 
goes to infinity at a finite altitude. On the other hand, 
if the vehicle orientation were such that .1 > C (like 
the vertical dumbell) then C,,, would be negative, and 
the period of the pitching mode would tend to the limit- 


r l B r 
= /3 4 ae C orbit 


which could be greater or less than 7 Gy nit. 

The phugoid periods obtained with the complete 
equations, cases (A) and (B), are the same—.e., the 
thrust law has no effect. It is interesting to note that 
this curve becomes asymptotic to the orbital time; 
above 400,000 ft they are essentially the same. This 
result is not observed, however, when the terms con- 
taining go are neglected. This approximation leads to 
an oscillatory mode up to 300,000 ft, but at the higher 
altitudes it degenerates into a pair of nonperiodic 
modes. The qualitative nature of the approximate 
solution is therefore wrong at the highest altitudes. 
On the other hand, the approximation 6 = 0, although 
it gives results which are wrong (by as much as a factor 
of 10) at intermediate altitudes, still has the correct 


ing value 


asymptotic behavior. 

The asymptotic behavior of the phugoid period may 
be regarded as a consequence of the equilibrium relation 
for the vertical force, the first of Eqs. (2.1). At low 
speed, we may neglect the term m*/ro, and the result 
is the classical, aerodynamically induced phugoid 
oscillation, in which the speed, flight-path angle and 
height vary in a characteristic manner. At the other 
extreme, in space where py) = 0, the aerodynamic lift Ly 
vanishes and the perturbed motion is that of a satellite 
in a circular orbit. The quantities speed, flight-path 
angle, and altitude then vary periodically with period 
equal to the time taken to describe the orbit, which is of 
course an ellipse of small eccentricity. The equations 
describing this case are derived in reference 8, and agree 
with the limiting results for the flight path in the phu- 
goid mode obtained from the full set of equations. The 
latter thus give the correct limiting solutions as the alti- 
tude tends to zero or infinity. 

Above about 400,000 ft, the X and Z aerodynamic 
forces are negligible, as evidenced by the fact that the 
phugoid period has already reached the empty-space 
limit. If we neglect the corresponding terms of Eq. 
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(2.17) and the damping terms Cx, and Cm,, which are 
also negligible, the characteristic equation of the set re- 


{* Cma Cm 
d (a: + JO. me ) =0 (4.2) 
ro~ lz 


The zero root is the limit of the spiral mode, the root 
+il/ro is the phugoid, giving a period equal to 7 ,,i:, 
and the last root is that of the short-period mode, cor- 
responding to Eq. (4.1). When the latter two roots are 


duces to 


equal—i.e , when 


(1/ro)? —(Cmg + Cig)/tg 


then the system is unstable, since the solution contains 
terms of the form / cos w/—t.e., a divergent oscillation. 
It may be expected that either large-amplitude pitching 
oscillations, limited by nonlinear effects, would develop, 
or possibly even tumbling. 


Damping of the Oscillatory Modes 

Fig. 9 presents the data on the damping of the two 
modes. There is only one curve for the pitching mode. 
All cases (A), (B), (C) and (D) gave the same result up 
to 400,000 ft. The damping of this mode is negligibly 
small at all heights in the range calculated, the number 
of cycles taken to reduce the oscillations to 1/2 ampli- 
tude increasing from 10 at 100,000 ft to about 5,000 at 
100,000 ft. This result has implications with regard to 
attitude control of such a vehicle: the contro] system 
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TABLE 1 

Altitude, Density, Uo, ee 

ft. X 1075 — slugs/cu. ft. m € Cw "a ft. /sec. 5 sec 
100 3.211 xX 1075 1.162 x 108 1.39 <X 1073 5.29 X 107 4.22 X 10. 3 5,940 0.230 22 , 550 
120 1.270 X 1075 2.935 X 103 3.53 X 10 5.74 <* 10? 2.76 X 1073 9,050 0.351 14,400 
160 2.4383 x< 107% 1.533 x 104 Loo x 7 8:70 x Ie7? 1.49 « 1073 16,800 0.652 7,860 
200 6.058 X 1077 6.16 xX 10° 7.44 x 10 1.00 <X 10= tie & 4°? 22 ,300 0.865 5,930 
240 1.268 <X 1077 2.965 X 105 3.58 X 107 7.66 X 107 1.002 <* 1073 24,900 0.967 5,320 
280 1.702 X Io- 2.19 XK 106 2.65 5.36 S.40 x. Or 25,600 0.995 5,210 
300 6.065 xX 107° 6.15 x 106 7.47 1.500 K 10 9.765 <x 10- 25,600 0.998 5,190 
350 4.957 XK 1071 7.52 X 10° 8.85 X 10 1.770 X 10° 9:.fa XxX 107° 25,700 1.00 5,180 
400 6.565 < 107?! 5.68 X& 108 6.66 X* 10? 1.332 X 10? 9.75 K 10 25,700 1.00 5,210 
450 titi xX 30-* 3.36 X 10° 3.925 X 10% 7.850 X 103 9.75 X 10 25,700 1.00 5,230 
500 2.862 * 10712 1.304 X 10% 1.522 X 104 3.044 x 104 0.75 xX lO 25,650 1.00 5,240 
600 4. 499 ¥ Jo-™ 8.29 X* 10" 9.66 xX 104 1.932 x< 105 9.82 X 10 25, 500 1.00 5,270 
700 ? x 107% 2.925 X 10" 3.39 X 105 6.78 xX 105 9.838 xX 10~ 5, 450 1.00 5,300 

itself would have to provide virtually all the damping. arrow), then a/@ = 0, |7/6| = 0.159, and # lags @ by 90°. 


In contrast to the pitching mode, the phugoid damp- 
ing is strongly affected by the thrust law, and by the 
approximations made in the equations. It is quite clear 
that the 8 and q terms must all be kept if correct results 
are to be obtained. The damping with an air-breathing 
engine is poor at all heights, whereas damping with a 
rocket engine remains good up to nearly 400,000 ft, at 
which height the period is already nearly equal to the 
orbital period. The large difference between the two 
thrust laws may be explained qualitatively as follows. 
The phugoid oscillation is damped mainly by the drag 
force, which varies with speed and height and extracts 
energy from the oscillation. When the thrust is con- 
stant, the net X force varies as the drag, and provides 
damping. However, when 7 and D both vary with p, 
then the effect of D is partially offset by the variation 
of 7 with p. Only the variation of D with speed now 
remains to provide damping. The heavy damping of 
case (A) at 300,000 ft appears to be a result of the peak 
in 6 at that altitude. 


Time Constant of the Nonoscillatory Mode 


Fig. 10 shows the characteristic time of the nonoscil- 
latory mode, with the orbital period plotted as well, for 
reference. This mode is stable for the air-breathing en- 
gine up to about 350,000 ft. Above this height it is un- 
stable, and values are not shown for this range since they 
exceed 10® sec. For the rocket engine, and for the 8 = 0 
approximation, the mode is unstable. When it is as- 
sumed that go = 0, this mode disappears and the charac- 
teristic equation reduces to a quartic. In the most un- 
stable case—a rocket engine at about 200,000 ft—the 
time during which a disturbance doubles is about 1/5 
of the orbital period—i.e., any disturbance will increase 
by a factor of about 2° = 32 in one revolution. 





The Mode Shapes (Eigenvectors) 

The characteristic relations among the variables in 
the periodic modes are shown in Figs. 11 and 12 (phugoid 
and pitching modes). The vectors drawn repre- 
sent complex numbers, the magnitude and argument of 
which correspond to the relative amplitudes and phase 
angles of the periodically varying quantities. 

When the longitudinal axis of the vehicle coincides 
with the flight path at all times (as in the ideal flight of an 


This is very nearly the case for the phugoid (Fig. 11) 

all three altitudes shown. Hence in the region of ap- 
preciable density, below 400,000 ft, the phugoid may 
be called a “trajectory” or “‘arrow’’ mode. The phase 
angle of 7 is seen to vary from nearly 180° at 100,000 ft. 
to 90° at 400,000 ft. In low-speed, low-altitude flight, 
0, the phase angle of 7 is close to 90°. The 
The dendiey 


or when 8 = 
departure found here is due to the @ effect. 
gradient provides in effect a spring stiffness with respect 
to height such that the associated perturbation in lift 
is always directed toward the reference altitude. This 
effect would by itself cause the maximum in # to occur 
when 7 is zero, or when @ is a minimum—i.e., the phase 
of # = 180°. On the other hand, the gravitational 
contribution to the phugoid tends to make the maximum 








Fic. 10. Time constant of nonoscillatory mode. 
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Fic. 11 Vector diagram of phugoid mode Rocket engine. 


velocity occur at the lowest point of the oscillatory 
path—i.e., when 7 is a minimum and @ is zero (the phase 
of u = 90°). The net result is a phase angle gener- 
ally between 90° and 180°. 

Above 400,000 ft, the two periods approach one 
another, and coupling effects are found in the attitude. 
Physically, the pitching moment C»,a + Cr is not 
strong enough to make the vehicle ‘‘track’’—+.e., it can- 
not produce the rotations necessary to suppress a. 
Fig. 13 is an example of the time history obtained when 
the two periods are not very different from one another. 
It shows the variation of a with time following an 
initially nonzero value of a, at an altitude of 485,000 ft. 
The classical phenomenon of beats is clearly shown. 
When the two periods are exactly equal, the system 
is dynamically unstable, and the simple picture of nor- 
mal modes breaks down. An example of this is given in 
Fig. 14, which shows an analog computer record of the 
oscillations which follow an initial disturbance in speed. 
The amplitudes of # and ? are constant, whereas those 
of a and @ are divergent. In the range of altitudes for 
which Eq. (4.2) applies, we obtain the following matrix 
equation for the mode shape: 


Dur —2¢ 0 0 ae] 
0 —Cmn, ipA 0 a/G6 
0 L/ro 0 nN q 6 
L/ro 0 l l/ro | 7/0 
iin 
= Cm (4.3) 
l/r 
A 
The solution of this equation for a/@, with 2e = C, 


(see Table 1), is 
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a 2? + 2(1/r%)? — (Cg lp) 


6 (Cm,/tp) + 2(l/ro)? 


(4.4) 


With the roots obtained from Eq. (4.2), Eq. (4.4) can 
be rewritten as 


= N= Ns — —- | (4.4a) 


6 (Cmg lp) — 2Xph = 


where A,.,. and A,». are respectively the short-period and 
phugoid roots. In the pitching mode, with A = A,.,., we 
get 

(a@/@)s.p l 


and in the phugoid mode 
Apa.” — Acy.” 


(Cma a = 2Ash.* 


(a@/) ph + | (4.5) 
At lower altitudes, —Cmg/tz > Cmg/tg, Azs.p.2 = Cmg/ip 
and Apn.” < As.p.”, and in that case (@/@),n. — 0, as in 
Fig. 11. In the neighborhood of the crossover altitude, 
however, when Aj». differs only a little from A,.,., sub- 
stantial values of a are developed in the phugoid motion, 
as shown in Fig. 14, and it can no longer be character- 
ized as an “‘arrow’’ mode. When the two periods are 
exactly equal, Eq. (4.5) yields (@/@)pn 1, which also 
is verified by Fig. 14. It is interesting to note that a 
possible condition is (@/@)pn. = ©-——1.e., 6 = 0. This 
occurs at the altitude where 


Cme 1B = 2X ph 2 


The significance of a small or zero value of 6 is that the 
vehicle z axis is always aligned with the local vertical, 
even though the flight path is an ellipse. This result 
might be useful in connection with attitude control. 

Fig. 12 confirms the preceding deduction about @,/@ in 
the pitching mode. It is seen to be close to unity even 
at low altitude. 








 « 
———————————— — (2) 
(a) ALTITUDE 100,000 ft. 
a, r negligible 
_ 0.x 





(b) ALTITUDée 300,000 ft. 
u, © negligible 
Rocket engine 


Fic. 12. Vector diagram of pitching mode 
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Fie. 13. Transient solution (linear) for @ following an initial 
condition on a. Altitude = 485,000 ft 


The facts that 7 and ¢ are negligible and 6 = a@ mean 
that the c.g. describes the undisturbed reference path 
at constant speed ms, while the longitudinal axis oscil- 
lates in pitch about it. This situation persists at all 
altitudes for which this mode is stable. 

The characteristics of the spiral mode are shown by 
the ratios of the variables, which for Case (A) are given 
in Table 2. 


TABLE 2 
Altitude, ft u/6 a/é 7 /0 
100,000 885 ma 87 
200, 000 168 =f) 1.86 
400 , O00 —32.8 ==) 67.6 


The vanishing of a in this mode, as in the phugoid, 
shows that the vehicle axis remains aligned with the 
flight path—i.e., this is an ‘‘arrow’’ mode. The flight 
path described is a spiral, proceeding away from the 
reference orbit. The ratio */@ can be shown, for a tra- 
jectory mode, to be simply related to the real charac- 
teristic root \ by 


#/6€ = L/ror 


and hence it simply indicates the rate of divergence, 





os | 
ne ee / \/ heiiienii 


Fic. 14. Transient solution following initial condition @ = 
—(.01—.e., 1 percent speed reduction. Altitude = 490,000 ft. 
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small 7/6 implying rapid divergence and vice versa. 
The variation of 7/0 with height is more interesting, in 
that it changes sign at some altitude /’ near 350,000 ft. 
This means, with the sign convention chosen, that in a 
descending spiral (7 and @ negative), the speed increases 
above h’ and decreases below it. This result is con- 
sistent with the behavior of an unpowered satellite 
spiralling in. The altitude 350,000 ft is also represen- 
tative of the height at which it ceases accelerating and 
begins to slow down. 


(5) Nonlinear Effects 


It was pointed out in Section 2 that the validity of the 
linearized solution at high altitude requires investiga- 
tion, since the condition g < q is easily violated. Con- 
sequently, a study has been made of the importance of 
the second-order terms. The terms gw and gu in the 
first two of Eqs. (2.8) were retained, and the set of equa- 
tions mechanized for analog computation. Fig. 15 
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Fic. 15. Example of second-order effects. Initial condition on 


a. Altitude = 100,000 ft. 


shows a record taken for 100,000 ft altitude and an 
initial value of a = 0.25 rad. The envelope of the 
second-order term ug is seen to be at most about 10 per- 
cent of a. Thus for this case, neglecting the second- 
order effect does not introduce a major error. This con- 
clusion was found to apply at all altitudes investigated. 
That is, for amplitudes of the disturbance quantities up 
to about 0.2, the second-order terms do not introduce 
important effects. The explanation for this will be 
given in a forthcoming UTIA report by Rangi.' 


Appendix: Alternative Description of the 
Dynamic Instability 


For the same conditions that lead to Eq. (4.2)—i.e., 
neglect of all aerodynamic terms except Cm,a—the 
following equation for the pitching degree of freedom 
(@) can be obtained without approximation from Eq. 
(2.17): 

(D? = New. 0 = (Cm ip)Y + (1 1) Dit =< vr) (A.1) 


where \,.,. has the same meaning as in Eq. (4.4), and y 


(Continued on page 832) 
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Satellite Motions About an Oblate Planet 


MAURICE L. ANTHONY* ann GEORGE E. FOSDICK** 
The Martin Company 


Summary 


The drag-free path of a vehicle in the vicinity of an oblate 
planet is discussed. Although the solution found is approximate 
in the sense that it is a truncated power series in the oblateness 
parameter J, there is no restriction on either the eccentricity or 
the inclination angle. Approximate analytic expressions are 
derived for the radial position, speed, angular momentum, devi 
ation from the initial plane of motion, and rate of apsidal ad- 
vance. In addition, the injection speed required for prescribed 
apses is derived. 

rhe differences between first-order results and inverse-square 
results are found from the analysis, and numerical results are 
given. A comparison of the analytical results with the results 
of integrating the equations of motion directly on a digital 
computer indicates excellent agreement over a wide range of 


inclination angles and eccentricities 


(1) Symbols 
Vo/V- 

D = coefficient of fourth harmonic in the gravitation po 
tential function for an oblate spheroid (approxi- 
mately 1.060 « 10~® for earth) 

r— 1% 


(J/6)( R? 


¥ << Te 


(J/6)( R?/re) 


G = Newton’s constant of universal gravitation 

grad = gradient of a scalar function 

i = orbital inclination to the equatorial plane 

J = coefficient of the second harmonic in the gravitational 


potential function for the oblate spheroid (approxi- 
mately 1.6388 « 10~* for earth) 


L = defined by Eq. (3.47) 

Ly = defined by Eq. (3.56) 

\f = mass of the oblate planet 

V = defined by Eq. (3.48) 

V4, = defined by Eq. (3.57) 

\ = defined by Eq. (5.4) 

eS = angular momentum per unit mass, to first order in 
J (= Po + JP) 

id = angular momentum per unit mass about a spherical 
earth (= Py = ro Vo) 

q = JR*/c'ro? 

= radial position of satellite 

ra = second prescribed apsidal radius 

ro = initial value of r (at an apse) 

re = radial position for a satellite of a non-oblate planet 

Vn =rtte-h=*f 

‘a = rat an additional apse 
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Zz, 


R = equatorial radius of the oblate spheroid (approxi 
mately 2.092601 & 10° ft = 3441.7779 nm for 
earth 

t = time 

u = 1/r( = wu + Ju, to first order in J 

U = gravitational potential function, Eq. (3.2) 

V = orbital speed 

V = V u r, circular speed at r = ”% for a non-oblate 
planet 

V = orbital speed for a non-oblate planet 

Vi = initial speed 

AVo = Vo — V- 

x) = inertial coordinate system with Y and Y in the equa 

Y? torial plane, and Z pointing to the north pole in a 

z\ right-handed system 

x] = intermediate coordinate system with x and y in the 

yp initial plane of motion, and z perpendicular to this 

| plane to form a right-handed system (x = X) 

a = angle between Z-axis and the radius vector 

€ = J(R/rn)? 

n = ¢c* — 1 (eccentricity, e = | 7 | for an elliptical orbit 

6 = angle between the z-axis and the radius vector, Fig. 1 
(= (r/2) + J6, to first order in J) 

6, = @for non-oblate planet = 2/2 

d = 6n/e = (2AV0/ V-)/(J/6)( R/1r)? 

im = GM, approximately 1.4076613 x 10'© ft.5/sec.? for 
earth 

é = independent variable for which the first-order analyti- 


cal results for rand V are periodic 


bo = initial value of 

He = £for whichr = r* 

p = rsin 6 

¢ = spherical coordinate (Fig. 1) representing the central 
angle in the initial orbital plane from the x (= X) 
axis to the projection of the radius vector in the xy 
plane (= — + J¢é, to first order in J) 

¢1 constant chosen to eliminate secular terms in r and V 

Gs = g fora non-oblate planet 

¢ = initial value of ¢ 


Aw = rate of apsidal advance, referred to inertial space 
A prime superscript indicates differentiation with respect to &. 
A dot over a variable indicates differentiation with respect to ¢ 


(2) Introduction 


Y  * seonse orem RESULTS are presented for the path of a 
satellite solely under the gravitational influence of 
an oblate planet. The path begins at an apse—e., 
where the radius and velocity vectors are orthogonal. 
The potential includes the second harmonic with the 
oblateness parameter J as a coefficient. An approxi- 
mate solution of the path equations is found following 
the method of Lindstedt!~* retaining only the terms of 
first power in J. No restriction is placed on either 
eccentricity or inclination angle. Thus, for example, 
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Fic 1. Coordinate systems and notations used. 


the analytical results may be applied to cases where 
the eccentricity is 0.5 (or to eccentricities exceeding 
1.0, in which case the vehicle is no longer a satellite*!®) 
for all inclination angles, including 90°. This paper 
constitutes an extension of several recent papers that 
treat motions with small eccentricities and often with 
small inclination angles, as exemplified by the work of 
Blitzer,>~? Weisfeld,> and Wheelon,®® Brouwer,’ 
Moulton,’ Orlov,'® Roberson,!! and King-Hele and 
Gilmore.'? For several different eccentricities and 
inclination angles, the paths of earth satellites have 
been found both by the approximate analytical results 
developed in this paper and by direct numerical integra- 
tion of the equations of motion. Comparisons made 
between the results of these two methods indicate 
that the analytical solution is correct to within 300 ft 
in radial position, and 0.1 ft/sec in speed. 

Section 3 develops the method for obtaining analyti- 
cal results pertaining to radial position, speed, angular 
momentum, deviation from the initial plane of motion, 
and apsidal advance of a satellite. It also includes a 
discussion of nearly circular orbits, the differences be- 
tween first-order oblateness results and inverse-square 
results, and the speed required at one apse to reach the 
second prescribed apsidal distance. Section 4 is a dis- 
cussion of these analytical results. Numerical results 
pertaining to these concepts are presented in Section 5. 


(3) Methods of Obtaining Analytical Results 


(a) Statement of the Problem 

The analysis presented in this paper concerns the 
determination of the path of a satellite under the sole 
influence of the gravitational field of an oblate planet. 
Initially, the position vector of the satellite is arbi- 
trary, and the velocity vector is orthogonal to the radius 
vector (so that the initial point is an apse) but is other- 
wise unrestricted in direction. Fig. 1 shows the initial 
plane of motion defined by the initial radius and ve- 
locity vectors [and denoted as the xy plane], where the 
x-axis is its intersection with the equatorial (X Y) 
plane of the planet. The Z and z axes are perpendicu- 
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lar to the equatorial (XY) and initial (xy) planes, 
respectively, so that both (X,Y,Z) and (x,y,z) form 
right-handed systems. The angle 7 between these 
axes is the inclination of the initial plane of motion 
with respect to the equatorial plane. Also shown in 
Fig. 1 are the spherical coordinates (7,6,¢) in terms of 
which the analysis is undertaken. 

If a planet were spherically symmetric, the force 
field would be Newtonian, and the motion would re- 
main in the initial plane of motion, determined by the 
initial radius and velocity vectors. For an oblate 
planet, the force field is generally neither central nor 
inverse-square nor is the path a plane curve. How- 
ever, for several revolutions the path remains very 
near the initial plane of motion (contrary to the path 
depicted in Fig. 1, which is exaggerated to indicate 
the angles). 

The equations of motion of a satellite in a potential 
field may be written in the vector form 


r = —grad U (3.1) 


For an oblate spheroid the potential U is given by 


Jeffreys'* as 


=H) 
U = R ( + J 3 cos* @ - + 


a " mie. soe 
— (35 cost a — 30 cos? a + 3) » (3.2) 
ma rj § 
where a is the angle between Z axis and the radius 
vector, and J and D are dimensionless constants that 
are empirically determined. For the earth these are 
given approximately by 


J = 1.638 X 10-4| 


D = 1.06 X 10-5 § (3.3) 


In the following analysis is found an approximate 
solution, correct to the first power in J; the higher 
powers of J are omitted because of their relatively 
small size. Since both D and J/* are of the second 
order in the planetary ellipticity, terms containing 
these factors are dropped from the approximate solu- 
tion. Such a solution is valid for several revolutions, 








Fic. 2. Symmetry of positions of Satellites 1 and 2 with respect to 
a vertical plane. 
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9 Pe — x 
Y = ycosi — gsini (3.5) 
Z=ysini+zcosi 
= NS > th I ial li 
—— a so tha » equatorial co ates may be expresse 
SS! —] (Sap=. a. | t the equatorial coordinates may be expressed 
| N LD in terms of 7, 8, and gas 
\ ~ ‘ s/f 
\ \ \ ne A 7 , . 
\ . \ 4 ff X =rsin@cos¢ 
\ X 7 if Y = r [sin @ sin ¢ cos i — cos @ sin 7] (3.6) 
\ " J, Z =r |sin @sin ¢ sin7 + cos 6 cos 7] 
x Phaee— f+ — The quantity cos a is given by 
K A cos a = sin @ sin ¢ sin 7 + cos 6 cos 1 (3.7) 
te Since in the (7,6,¢) system the gradient of a scalar 


. has the components 
1G 


Effects of oblateness on orbital radius for several values 


of & (7 = 0.1, 7 = 45° ol’ 10U’ 1 OU 
grad U = , oe (3.8) 
_Or +r OO rsin6@0¢ 


until the first-order term in J no longer dominates the 


; ane > components of acceleration are 
perturbation. 1 the components of acceleration are 


The transformation relating the (x,y,z) system to the - g i a 
(r,0,0) system is r= ty — 0 — re" amr €, rw) = 
, i r dt 
x =rsin@cos¢ —_—" 
y=rsin@sin¢ (3.4) rg” sin 6 cos 6, — (r°¢ sin? @) (3.9) 
: rsin 6 dt 


z= rcos#@ 
and that relating the (X,¥,Z) system to the (x,y,z) ™ teeta — y negpirod bog ne 
combining Eqs. (3.1) and (3.2) with D = 0, and Eqs. 


system is psaggse 
(3.7), (3.8) and (3.9): 


r — 16? — rg? sin?@ = —(u/r?) {1 + J(R/r)? (1 — 3 cos?a)} 
(d/dt)(r?6) — (1/2)r?¢? sin 20 = —(uJR?/r*){ (sin?i sin’g — cos*/) sin 26 + sin 27 cos 26 sin ¢} (3.10 
(d/dt)(r2¢ sin’@) = —(uJR2/r*)} sin’ sin? @ sin 2¢ + (1/2) sin 27 sin 26 cos ¢} 


Although in this paper we present the solution of the general problem for an arbitrary initial position, to facilitate 
the discussion of the method of solution we adopt here the initial conditions that at t = 0 the satellite is in the equa- 
torial plane, moving horizontally with speed V)in the xy plane. Thatis, 

r=n, 0=7/2, ¢=0 *+=0, 6=0, ¢ = Vo/to (3.11) 

The solution of Eqs. (3.10) subject to the initial conditions (3.11) determines the manner in which the polar coordi- 
nates (r,0,¢) depend on time. The equatorial coordinates (X,Y,Z) may then be determined as functions of time by 
Eqs. (3.6). However, even in the classical Newtonian case (two-body problem) where J = 0, it is difficult to deter- 
mine the coordinates as functions of time and is not as informative as it is to find the equations of the path. Asin the 
classical treatment of the two-body problem, the time ¢ is eliminated to obtain the equations pertinent to determining 
the path of the satellite. The time may be recovered at a later stage and expressed in terms of the position of the 
satellite, in a manner similar to the treatment of the two-body problem. 

(b) Determination of the Path and Speed of a Satellite 


The path of motion is determined by the introduction of 


u=i1/r (3.12) 


P=re = ae (3.13) 
Thus, for example, if Sis a function of g, which is in turn a function of ¢, 
dS/dt = (d¢e/dt)(dS/d¢e) = u?PdS/de¢ (3.14) 
By eliminating all derivatives with respect to time in this way, Eqs. (3.10) may be written 


(d/dy)(P du/dg) + Pu[(d0/dg)? + sin’@] = (u/P)}1 + JR2u*(1 — 3 cosa) } 
(d/dy)(P dé/dg) — (P/2)sin 20 = —(uJR*u/P)} (sind sin’¢ — cos’) sin 26 + sin 27 cos 26 sin ¢} (3.15) 
(d/dg)(P sin) = —(uJR2u/P)}sin2i sin’@ sin 2¢ + (1/2)sin 27 sin 20 cos ¢} 
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The motion remains in the plane established by the 
initial radius and velocity vectors if and only if @ = m 2 
throughout the motion. Direct substitution into the 
second of Eqs. (3.15) reveals that @ = 7m, 2 is a solution 


if and only if J sin 27 = 0. Thus, plane motion can 
occur in three distinct ways: 
(1) J = 0, and 7 is arbitrary —the classical problem 


of the inverse-square central force field solved by 
Newton. 

(2) ¢ = 0, and J is arbitrary —the problem of motion 
in the equatorial plane. The force field is still central 
because of the assumed symmetry in the figure of 
the planet, but it is not inverse-square. 

(3) ¢ = w/2, and J is arbitrary the problem of 
motion in the polar plane. The force field is neither 
central nor inverse-square. 

The first case, which is free of oblateness effects, 
can be solved in closed form and can form a quali- 
tative frame of reference for the more general 
problem. The second case, equatorial plane motion, 
is of interest both in connection with establishing an 
orbit in which maximum use is made of the rotation of a 
spheroid about its axis, and in connection with a 
satellite that is constantly directly over a fixed point 
of the spheroid—-e.g., a 24-hr. satellite of the earth. 
The third case, polar plane motion, is of interest in 
connection with the Discoverer series of satellites. 
In addition, the two cases of plane motion, : = 0 and 
1 = 7/2, are intimately connected with the problem of 
general motion in the vicinity of the spheroid, as is 
shown in Section (4a). 

An approximate solution of the differential equations 
(3.15) is found by the method of Lindstedt!~*—1.e., 
by assuming power series expansions in the small param- 
eter J for all variables and then truncating the series 
beyond the first power of /. Thus, the differential 
equations are not satisfied exactly. However, substi- 
tution of the approximate solution into the equations 
of the path leads to expressions of the form 


(zero) + J(zero) + 


J*(a function of variables and J) = 0 (3.16) 


Since J is a quantity small compared to one, the differ- 
ential equations are “‘nearly’’ satisfied. 
the solution is approximate. The truncated power- 
series expansions referred to may be expressed as 


u = Uu(E) + Jm(é) 
P — P(&) -t- JP \(&) tachie? 
6 = (r/2) + J6,(E) 


where the new independent variable & is defined by 
¢g¢ = &(1 + J¢1) (3.18) 


where ¢; is a constant that will be chosen later so as to 
eliminate secular terms in 4 expressed as functions of &. 
The quantities 1) and /y are the inverse radius and the 
angular momentum per unit mass for the two-body 
problem, in which J = 0. The quantities 1, P:, 6, 
and ¢; all have the small multiplier J. Hence u, ?, and 


In this sense 
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Fic. 4 Effects of oblateness on orbital radius for several values 


of & (yn = 0,1 = 45 


§ are predominantly given by the solution of the two- 
body problem. Such small perturbations as Ji, JP, 
and J6,; are the concern of this analysis. 

The length of the projection of the radius vector on 
the initial plane of motion, p, is the same as r, to the 
first order in /, as can be seen from 


p =rsin@ = rsin[(x 2) + J6| 


r+ O(J*) (3.19)* 
Aiso, the angular momentum per unit mass is /?, to the 


first order in J: 


pe = re sin’?é = rg + O(J?) = P+ O(J*) — (3.20) 
Thus, 7 and P may be thought of as the radial coordi- 
nate and the angular momentum, respectively, of the 
particle, as if it were moving entirely in the plane 
established by the initial motion. This conclusion 
may also be reached by observing that the first and 
third of Eqs. (3.15) are not affected, to the first order 
in J, whether 6 = 7/2 or 6 = (7/2) + J0\(é). These 
two equations, therefore, are only weakly coupled 
with the second of Eqs. (3.15). They can be solved 
as if the entire motion occurred in the initial plane; the 
second of Eqs. (3.15) may then be solved to determine 
how much the motion departs from this plane. This 
whole process is correct to the first order in J. 
Mathematically, the solution is obtained by substi- 
tuting Eqs. (3.17) and (3.18) into the equations of 
motion (3.15), expanding all trigonometric functions 
in powers of J, and dropping terms higher than the 
first power of J. These operations result in three 
equations for the determination of the five unknown 
functions (i, %, Po, Pi, 6:) and the constant ¢;. De- 
noting differentiation with respect to & by a prime, 
the third of Eqs. (3.15) becomes the following identity 


in J: 
Rar + JP,’ = giPo’ + 


[(ukR? sin?2)/Poluy sin 2 = 0 (3.21) 


* Terms of second order in J are indicated by O(/? 
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Since this relation should hold independently of the 
value of J, the coefficient of each power of / must be 


zero; that is, 
fy = @ (3.22) 


>] 


(uk? sin?7) Pol sin 2E = 0 (3.23) 


P,’ 2 


Simularly, the first and second of Eqs. (3.15) imply 


ty” + to = p/Po? (3.24) 
” : MY 2P; 
+ iy —24tty + - 1291 — 
P,.? ( P, 
R?[uy? + sin? @(uoty’ sin 2E — 3x9? sin? &) r 3.25) 
0” + 6 (ulk?uy Py”) sin 27 sin £ (3.26) 


Eqs. (3.22) through (3.26) are sufficient to determine 
the functions u, ?, and @ to first order in / provided 
initial conditions are appended. The initial conditions, 
given in Eqs. (3.11), may be expressed as follows: at 
f= () 


‘ ’ 


g 0 


“uy = P, le’ = tty’ 6; = 6,’ 04 


The solutions of Eqs. (3.22), (3.23), and (3.24) sub- 
ject to the initial conditions (3.27) are: 

P5. = Vo (3.28) 

sin77[1 + (4/3)n — neosé — 


P, — (1% Vo/2c4) (RK, 19)? 


and 


. 9) 
2 3 


n* l | ok 2 l pare i a per 
B (1 + ) + (- —-- yy > ”') sin*? + > (9° + sin*t) cos 2 + -n sin*7 cos 38& + -n° sin*2 cos 4€ (3.5 
Z 7 3 S 


All the terms in B give rise to particular solutions that either are constant or are periodic in &. 
ticular solution arising from .1 would contain the secular term & sin &. 
In particular, 7 = 0 will eliminate secular terms. 


will consequently make , a periodic function of &. 


TE MOTIONS 793 


_e ry 
‘Sa ee a 
Xo: ~ * 
YO — 
Waee 
. — 
N= 
hic. 5. Effects of oblateness on orbital radius for several inclina- 
tions (yn = 0.1; & = 0°, 905 
cos rE —m te $3) cos SE (3.29) 
TT (1/roec*)(1 + cos & 3.30) 
where 
y= c*— 1 = (V,/V.)? — 1 (roVo?/u) — 1 (3.31) 


) 


The use of these results in Eq. (3.25) leads to a differ- 
ential equation of the form 


where 


A = nl2 — 3 sin? — 2eict*(r, R)*| cost (3.33) 


) ~ 


However, the par- 
Any choice of parameters that makes 4 = 0 
Eliminating 


secular terms may also be possible by choice of / after yg; is chosen; however, these are very special cases that put 


some restriction on the problem. 
putting any restrictions on the parameters 9, ry, X, or 1. 


1 


With the secular terms in & thus eliminated, the solution for m4, subject to the initial conditions (3.27 


¢ 


= (1 ¢4)(R ro)*{1 — (3/2)sin*2 (3.35) 


On the other hand, the choice of ¢g; may always be made so that 4 = 0, without 


This choice is 


» 1S 


l R 2 4 7 l } on iz - 2 9 2 ; ; 
uy co Ny ) > > + ie ae ae sin*z | + | — I+ 3” + eats? sin*z | X 


rot 


cos £ — — (n* + sin*7) cos 2 — 7 n sin*7 cos 3& — ry n°” sin*? cos 4E 
) 


(3 


The solution of Eq. (3.26) subject to initial conditions 


(3.27) is 


sin 27 /R\? os : 
i, = 1(3 + 2n) sin & — 


6c4 r 
n sin 2 — 3& cos —{ (3.37) 


In spherical coordinates, the speed of the satellite is 


5 ] 
(3.36) 


given by 
V2 = r2 + (rb)? + (r¢ sin 6)? 
= P2\y?[sin? @ + (d0/dg)*| + (du/dg)*} (3.38) 
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: i ee. SO TABLE 1 _ = 7 
: oa ae General First-Order Results (& = 0 
2 
- 
= jl oe (2 - 3 sin” i) F (3.40) 
oe ir 
Ze Oo 
—— 2 
6 = is = ri £.. (3 + 2n) sin & - n sin 2€ - 3& cos '3.41) 
2 4 i \ | s c 
6c oO 
2 
2 4 , 1 J R . . 
P=rfr ¢ - - * 4 r Sin FA [3 + 4n - 3y cos ¢ - 3 cos 2¢ es 
6c oO 
- qj cos a 3.42 | 
2 
se lL qcsce+— oe L 13.43] 
r ee = 4\r . . 
rc 24c oO 
fe) 
l+1 wt he L ee 
= = l 1 - ‘OS ft 7” 24 rr 2 13.44) 
' ee ” Oo (1 + yn) (1 + 4 cos £) 
2 
V 2 
2 J : 
vf «ak. +n + 2n cos —& + a 7 [3.45 
4 4\r . ’ 
e 24c fe) 
2 - 
* y1 + + 27 cos & — rR ¥ M me 
V= ex -* nrg 5 [3.46] 
: ) (1+) (1+ + 2n cos €) 
where 
n..* . = oe 7 
L =] 24 + 12H + [- 12 + 32y - 15y sin” i] + |- (24 + 8y + 
2 Z Z Zz 
+ (16 - 2747 + 16n ) sin’ i] cos € - 4 (1 + sin i] cos 2& - 
2 Z 
- 55 sin’ i cos 3 - 7 sin i cos 4 [3.47] 
; 3 —— . e P 
M = 16 3- 3H - HRP +L - 3 + 4y - 3H sin i] + 2 2(- 12 + 12n - 
2 3 2 3 2 
- 4n + 3y ) + (16 - 45y + lon - 3y Sin i] cos & + 
2 2 3 3 2 
+ 16 (+ + sin i) cos 2E + [4n + ( 26 + 3y ) Sin i} cos 3t + 
Z Z 3 2 
+ 16 sin i cos 4& + 3y sin i cos 5¢ (3.48) 











By using 0 = (7/2) + J6:, and introducing é, the fol- 
lowing expression, correct to the first order in J, is 
obtained: 


V2 = Pn? + (1 — 23 ¢1)(u’)?] (3.39) 


Alternatively, V* may be determined from conserva- 


tion of energy. This result and others for equatorial 


= By use of the 


launch (&) 0) are shown in Table 1. 
above methods, analytical results may also be obtained 
for the path of a satellite that is initially at an apse 
not located in the equatorial plane (£5 # 0). These 
2, in terms of (£ — 


a) 


results are presented in Table £). 
The initial value of — (denoted by &) may be obtained 
from Eq. (3.49) when the initial central angle, ¢o, is 


specified. 
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For a nearly circular orbit—1.e., an orbit for which 7 
is of the order of J—the omission of products of order /J* 
leads to the results of Table 3 (for & 0) and Table 4 
(for &, # 0). 


(4) Discussion of the Analytical Results 


Several interesting features of the solutions just pre 
sented are discussed in this section. Included are 
(1) a few remarks concerning the deviation of the path 
from the initial plane of motion, the rate of epsidal 
advance, the size of the parameter (J? c'r)*), a super- 
position principle, and some properties of the solution 
that follow from the symmetry of the force field; (2) 


a discussion of the positions at which epses occur, which 
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shows that a circular path may occur only in the equa 
torial plane, that in general a satellite path is not sym 
metric with respect to its apses, and that two, three, or 
four apses may occur in a single revolution; (3) a de 
termination of the initial speed required at one apse to 
reach a second prescribed apsidal distance. 


(a) Properties of the Solution 


The result given in Eq. (3.50) shows that the motion 


remains plane if 7 0 or r/2—1.e., for equatorial and 
For other angles of inclination there is a 
secular growth of @ with (é — &,); 1n addition, this devi 
ation from plane motion is maximum for 7 


The functions 7, P, 7, and |'all have the period 27 in 


polar planes. 


w/4. 


the varizble ( — &,); Eq. (3.49) implies that they have 


TABLE 3 


Nearly Circular Orbits (& 











| 
[ J/R \ 2 
= ‘ + 3 (r (2 - 3 sin i)] & [3.58] | 
O | 
ce J oR - 
2= + > te sin 2i (sin — - & cos &) [3.59] 
i a ae 
2 
tit = - 
|; P= 7. 1 - = ( - sin i (l - cos 2€) [3.60] 
= O 
i 2 | 
u=- L=- {i = cos §€) + Ji = [ca - cos €) + 
ma O 
==. 
+ % sin i (- 3+ 4 cos & - cos 22)| [3.61] | 
R 2 
eee 1+ n(l - cos &) - J ay fa - cos —&) + 
O 
1 2 
+ vA sin i (- 3 + 4 cos & - cos 28) [3.62] 
2 2 . 
Vo= 7. 1 - 2n(l - cos —) + 2J ~ jc - cos —) + 
O 
aw 
+ 3 «sin i (- 3+ 2 cos —€ + cos 2€) [3.63] 
; R \" 
|e Vo Ll - n(l - cos £) + Ji — fc - cos —&) + 
| O 
i 2 : 
| + 3 sin i (- 3 + 2 cos —§ + cos 28) [3.64] 
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the following period in the space angle (¢ — go): 
T = 2x[1 + (J/2c*)(R/r)?(2 — 3 sin%z)] = 
249 + Aw (4.1) 


Thus, for example, the apses will advance in inertial 
space by 


Aw = (rJ/c*)(R/1)?(2 — 3 sin*z) rad./rev. (4.2) 


The maximum advance occurs when the satellite moves 
in the equatorial plane, in which case 

Awmar = (29J/c*)(R/r)? rad./rev. (4.3) 
The maximum regression of the apses occurs when the 
satellite moves in the polar plane, in which case 


Awmin = —(mJ/c*)(R/r)? rad./rev. (4.4) 


which is numerically one half the advance in the equa- 
torial plane. There is one inclination angle for which 
no advance of the apses occurs in inertial space. This 
angle is given by (2 — 3 sin*z = Q), the solution of 
which is approximately ¢ = 54.735°. 
A parameter of interest that appears in all the results 
is 
q = JR?*/c'ro (4.5) 


Because c may be small—e.g., for a path beginning at a 
large apo-apsidal distance and having a small peri- 
apsidal distance—for some paths, it appears that g may 
be large. However, by the use of Eq. (3.53), it can be 
shown that c‘ro” is larger than (1 + In |) *%min2- Since 
for a satellite path 7,,;, > R, the parameter gq satisfies 
the inequality 
q<J/A+ |n|)?<J (4.6) 

This relation is valid for satellites regardless of the size 
of ¢. 

The variables 9, P, r, and V (but not @), given in 
Tables 1, 2, 3, and 4, are of the form 


yi) = F+G sin (4.7) 


where F and G are functions of the parameters 7, 70, 
&, RX, J and the independent variable ( — &). Thus for 
any fixed set of these five parameters and (& — &), Eq. 
(4.7) implies that the following superposition holds: 


y(t) = y(0) cos*t + y(r/2)sin*z (4.8) 


Thus any of the quantities (radius, angular momen- 
tum, speed, or space angle) may be determined by 
forming the above linear combination of these same 
quantities which are solutions for the equatorial plane 
and the polar plane, using the same values of 7, fo, £6, 
R, J, and (€ — &) in each of the quantities. Thus, 
the motion of a satellite initially inclined to the equator 
is intimately connected with the motion in the equa- 
torial and polar planes. 

In general, the functions 7, lV’, and P depend on &» in 
two ways: (1) through the independent variable ex- 
pressed as (£ — &); and (2) through the coefficients 
that contain sin 2é and cos 2. However, for motions 
in the equatorial plane (¢ = 0), the coefficients no 
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longer depend on &, so that 7, V, and P depend on &) 


only in the combination (€ — &), as is readily seen 
from the rotational symmetry of the assumed force 
field. 


Because the coefficients in 7, V, and P contain & 
only through sin 2 and cos 2£), they are periodic func 
tions of & with a period of 7. Hence, the functions 
r, V, and P are the same functions of (¢ — &) for two 
satellites that have identical initial conditions except 
that one begins at & = 8 and the other at é& = mw + 8, 
where @ is arbitrary. This result may be expressed as 


r(éo, & = &) = 1(7 + £, & ae &o) (4.9) 


where the first entry indicates the argument of the 
coefficients, and the second indicates the independent 
(running) variable. Since 7, V, and P are periodic 
functions of ( — £) with period 27, the entire behavior 
of these quantities may be studied by examining the 
range 0 < (€ — &) < 2m. Eq. (4.9) is useful because 
the influence of the position (&) of the initial apse may 
be determined completely by considering only the 
range 0 < & < zm. 

A second useful relationship is due to the symmetry 
of the force field with respect to a plane through the 
polar axis of the spheroid. Fig. 2 shows two satellites 
in the same initial plane of motion. Since the force 
field has mirror symmetry with respect to the vertical 
plane, if the initial conditions also have this property, 
the satellites will always maintain a mirror image 
relationship. Mathematically, this relationship may 
be expressed as 


rl£o, E — £] i r(x = £0, 2r — (é = £) | (4.10) 


which implies that only the range 0 < & < 7/2 need 
be studied in order to determine completely the influ- 
ence of the position (&) of the initial apse on the radius. 


(b) Positions at Which Apses Occur 

Included in this section are discussions which estab- 
lish the following conclusions: (1) a circular path is 
possible only in the equatorial plane; (2) an apse occurs 
at (£ — &) = m only if the initial apse occurs either in 
the equatorial plane (& = (0) or at maximum latitude 
(&) = w/2), and that the functions 7, V, and P are even 
functions with respect to these two apses; (3) for speeds 
near circular speed for an oblate planet, the path may 
have two, three, or four apses per revolution; (4) for 
speeds not near circular speed, the angular position of 
the second apse is determined as a function of &, 7, ro, 
and 7. 

An apse is a point on a trajectory for which * = 0. 
The radius at which an apse occurs is called an apsida! 
distance, and the angle subtended between two suc- 
cessive apses is called an apsidal angle. If a particle 
moves in a circle, then every point is an apse, all apsidal 
distances are equal, and the concept of apsidal angle 
has no meaning. In cases other than this special case, 
if a particle is in a central force field in which the force 
intensity is a single-valued function of the radius, the 
path has at most two apsidal distances and one apsidal 
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angle. This well-known result for central force fields 
follows from symmetry of the path with respect to the 
apses, and is applicable to equatorial motion. In 
general, however, apses may occur with or without 
symmetry for a noncentral-force field, as in the case 
of nonequatorial motion. To determine where apses 
occur in this problem, the zeros of u’ are found. For 
all trajectories for which the angular momentum P is 
not zero, u’ = 0 is equivalent to * = 0. 
From Eq. (3.52) the derivative of u is found to be: 


du l j in (é bo) 
= + «aS ee 0) 
dé roc? \ . ; 
/ (4) ou in w(€ — &) + 
Ag Sin NE — &o) 
24c4 \roJ n=1 
B,, sin 2& sin®z cos n(& — &) It (4.11) 
where 


A, = 24 + 8? + sin*7[(20 + 


27n — 4n”) cos 2& — 36 — 12n°] 


Ay = 4)2n? + 
<b , — (4.12) 
sin’z [(2 + 3m?) cos 2% — 3n?]} 
Az = 15m sin*z cos 2& 
A, = 4n? sin?2 cos 2& 
B, = —(8 + L5n + 16”) B; = 15) (4.13 


By = 4(2 + 3n”) B, = 4n? 4 


If the oblateness parameter J were zero, as in the 
inverse-square force field, the path would be circular 
and hence “’ would always be zero provided n = 0. 
For any other initial speed (y # 0), apses occur at 
solutions of sin (E — &) = O0—i.e., (& — &) = Aaw—where 
K = 0, 1, 2, etc. These basic results for the inverse- 
square field evoke the following questions concerning 
satellite paths near an oblate planet: (1) under what 
conditions, if any, is a circular path possible? and (2) 
under what conditions are (£ — &) = Az apses of the 
path? 

(1) Circular Path—In order to have a circular path, 
u’ must be identically zero. Since g < J, no circular 
path is possible unless 7 is of order J. For if » were 
large compared to J, then the term —7 sin ( — &)) 
could not be cancelled by the remaining terms, which 
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are multiplied by J. As is consistent with neglecting 
terms of order J*, terms containing Jy, n* and higher 
powers will also be omitted in this discussion. The 
expression for u’ is therefore simplified to: 


du ere Maes : 
= — {Ci sin (€ — &) + Cosin 2(& — &) + 
dé 3ro 
C;[cos 2(& — &)) — cos (E — &)}} (4.14) 
where 
GC =3 "7 + « [ + sin*7 X 


( cos 2& — = (4.15) 
6 2/ \§ 


C. = € sin*i cos 2& 
C; = €sin*z sin 2% 

Since the three trigonometric functions of Eq. (4.14) 
are linearly independent, u’ can vanish identically 
only when C; = C, = C; = 0. The solutions of these 
three equations require that sin 7 = 0 and n = e€; & 
remains arbitrary. Thus a circular path is possible 
only in the equatorial plane, in which case the required 
speed may be found from n = e. 

(2) Apse at ( — &) = m—For an apse to occur when 
(¢ — &) = a, the derivative u’ must be zero. Direct 
substitution of (£ — &) = minto Eq. (4.11) yields 


du 7 2] (< 


= ; (1 2n?) sin’s sin 2¢ (4.16) 
dé 3roc® ) = ” f 


To 


Thus, in general, an apse does not occur at (& — &) = 


a. However, one can occur at (§ — &) = wifi = 
0, = 0, or & = 2/2. 
(a) Case 1 (i = 0): For motion in the equatorial 


plane, Eq. (4.11) may be written 


du sin (& — &) f 


dé = roc? = » t 
J {R\? e : { > 
[(3 + 7) + 2n? cos (E — &)] f (4.17) 
3c! ro 


Apses occur at (& — &) = Kw because the factor 
sin (£ — &) = 0. Additional apses may occur if the ex- 
pression in the brace is zero. Again, since gq < J, 
7 must be of order J to have this occur, in which case 
products of J and 7 are negligible. The brace vanishes 
only for » = ¢, which is the condition for a circular 
path (previously discussed). 

(b) Case 2 (& = 0): In this case, initially the 
satellite, but not the velocity vector, is in the equa- 
torial plane. Eq. (4.11) now becomes 


du sin & | J {R\? 
” | + ° x 


dé roc? : 


[1 — (2/3)sin*2(1 — cos &) + O(n) If (4.18) 


From the first factor, apses occur at £ = Kz, where 
K = 0, 1, 2, etc. The second factor may yield addi- 


-_ id } 


tional apses if the parameters satisfy the relation 
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n = e[1 — (2/3)sin? 7(1 — cos é*) | (4.19) 


where () < £* < a is the position where an additional 
apse occurs. A second additional apse will occur be- 
tween mw and 27. For & (), the functions 7, V, and 
P, given in Table 1, are symmetric with respect to the 
apses at both € = Oand & = r. 

(c) Case 3 (& = 2/2): In this case, the initial 
velocity vector is parallel to the equatorial plane. 
As in the previous case, apses occur at (£ — &) = Ara, 
where K = 0, 1, 2, ete., and additional apses may occur 
provided 

n = e}1 — (1/3)sin’/[7 + 2 cos (E* — &)]f (4.20) 


Furthermore, the functions 7, V, and P are even func- 
tions with respect to the apses at (& — £) 0 and at 


(— — &) = m. 
(3) Apse at (& — &) # mw—The previous discussion 
showed that an apse does not occur at (& — &) = @ 


except for certain initial conditions. The apsidal angle 
in the general case will now be considered. The dis- 
cussion consists of two parts: (1) paths for which 
n is large compared to J, and (2) paths for which 7 and 
J are of the same order. 

(a) Case 1 (\n| > J): In this case q is not of the 
order of J, so that additional apses are not encountered, 
and the path should be quite similar to the elliptical 
path of the inverse-square field. To determine the 
approximate location of the second apse, the following 


truncated power series is assumed: 


E—-f=ar+JA (4.21 


Ss 


where 4 is determined so that uw’ is zero to the first 
order in J. Substituting this value of ( — &)) into Eq. 
(4.11), setting uv’ = 0, and retaining terms to first 
order in J yields a linear equation in .1, the solution of 
which implies that the second apse occurs at 
(£ — &&) = wr — (2/3n)(JR2/c4ry?) X 

(1 + 2n7) sin*z sin 2& (4.22) 
This result reiterates the fact that an apse occurs at 
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Fic. 6. Effects of oblateness on orbital radius for several inclina- 
tions (yn = 0; & = 0°, 90°). 
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y 
£ y 
4f A - 
— ~ ee 
Nic. 7. Radius variations in nearly circular orbits for several 
initial speeds (2 = 90°; & = 0°, 90°). 
(E — &) eit = 0, & = 0, or & = 2/2. It also 


shows that for the range 0 < & < m2, the apsidal angle 
is smaller or larger than mz according to whether the 
initial apse is a peri-apsis or an apo-apsis—i.e., 7 > 0 or 
n < 0. The absolute value of the deviation of the 
apsidal angle from a is a maximum for the case 7 


1/2, & = m/4. Use of Eq. (4.6) in this case results in 
i 2J (1 + 29? 2] 

= se © eee =e (4.25 
Bln! + nl)? > 3ly 


Since in this case || is large compared to J, the devi- 
ation of the apsidal angle from 7 cannot be large. 
However, as |n| becomes sinall the deviation of the 
apsidal angle from m increases. 

(b) Case 2 [yn = O(J)]: For a nearly circular initial 
speed, the path may have additional apses. No apse 
can occur near (£ — &) = mw. Eq. (4.14) applies, since 
terms of order Jn and n° are negligible in this case. 
Examination of (4.14) reveals (1) that apses always 


occur at ( — &) = 2Kz, where K = (0), 1, 2, etc.; (2 
that if sin’ sin 2é = 0, then (§ — &) = w + 2K7z7 are 
also apses; and (3) that if sin’? sin 2& #4 0, then apses 
occur at (€ — &) = 2A and at any point where 

| jo q) 

aoa = 49127 < wie DE a 4 

(n/€) l = sin’ ba Oo 2& pe 

> -” -_ 

2 cos 2& cos (E* — &) + sin 2& X 
cos 2(¢* — fo) — cos (&*— &) TL yg, 
sin (&* — & f ins 


Thus, for any given position and direction, the speed 
may be determined by Eq. (4.24) so that an apse occurs 
at any preselected angle [except near (&* — &) = 7}. 
Depending on the value of » found in this way, 
there may be zero, one, or two more values of (* — &) 
that satisfy Eq. (4.24) and hence define more apsidal 


positions (see Fig. 9). 


(c) Injection Speed Required for Specified Apsidal 
Distances 


Given the location of an apse (é) and the direction 
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of the velocity vector, the path is then determined by 
the speed. Similarly, the initial speed is determined 
by the initial data r), 7, & and the prescribed second 


3(1 + )* \ro 


where all the parameters are known except 7. 


power series in J. The result is 


ry — to 2S R\* 9 ahd 
n 7: —?r,- sin-l cos 2& 
YT, t+ 1% 37 4(%4 + 70) \" ' 
or equivalently, 
: 27 4 J (R\*s i ; 
c= — Se eee 4 — 2ro? sin7 cos Zé, 
if “+ To ars ro 


As is the case for all the results found in this paper, for / 


(5) Typical Numerical Results 


This section contains the following numerical results: 
(1) effects of oblateness on orbital radius as functions 
of 1, &, initial speed, and (& — &); (2) relations between 
n, €, 1, and & to assure an apse at (€* — &); (3) varia- 
tions in 6as a function of & for & = 0 and several values 
of n; (4) the injection speed required at perigee for a 
specified apogee. 

Two functions pertaining to radial position are 
introduced 


PF 


(J 6)(RR? 7) 
r—n 


and see 
(J/6)(R2/ro) 


where r is the radial position derived in this paper, and 
r, is the radial position that would occur under the same 
initial conditions if the force field were inverse-square. 
For nearly circular orbits the function / focuses atten- 
tion on the small perturbations induced by oblate- 
ness, whereas these would be masked by the large 
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Fic. 8. Radius variations in nearly circular orbits for severa 
= , = ()°, 90° 


initial speeds (7 = 45°; & 
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apsidal radius. 


The speed is found from Eq. (3.52) where (& — &) is 
given by Eq. (4.22). If the second prescribed apsidal 


distance is denoted by r,, Eq. (3.52) becomes, to first 


order in J, 


sin*z [(5 + Sn — n*) cos 2é 3(3 + we y]t (4.2 


This relation is solved for y to first order in / by developing as a 


"> “J 
(74° + rato + P07) , sane l cos 2&5) ( 1.20) 
*) ae 
» £ , y ,{ - <9 4 Je { 7 
(v74° t+ Lalu T+ To i sin*7( 1 cos 2£) f (4.2% 
») 


() these results reduce to the inverse-square field results 


changes in r, if the orbit were eccentric. Hence the 
difference between the oblate solution and the inverse- 
square solution is utilized in presentation of numerical 
results regarding orbital radius. The quantity /; is 
identical to f when 7 = (0. Extensive graphical cover 
age of equatorial and polar orbits is given in previous 
reports!’ © for both orbital radius and speed. 

Numerous checks have been made between analytical 
results of the first-order theory and results obtained 
from direct numerical integration of the equations of 
motion. The agreement between the results of these 
two methods is excellent. One check was for the fol 
lowing case: 

» = 0.5, l 60°, £y () 


Initial altitude 300 st.mi. 


The results are shown in Table 5. 

The effects of oblateness on orbital radius, /;, as a 
function of (£ — &) are presented for (1) several values 
0 (Fig. 


of & at 7 45° for 0.1 (Fig. 3) and 9 
4+): (2) several inclinations at & 0°, 90° for 7 0.1 


(Fig. 5) and 0) (Fig. 6). Also, radius variations 


Fic. 9. Relations between 7, €, 7, and & to assure an apse at a 
specified (&* — & 








802 JOURNAL OF THE AEROSPACE 


in nearly circular orbits, /, for several initial speeds 
at & = 0°, 90° are presented in Fig. 7 for 7 = 90°, and 
Fig. 8 for 1 = 45°, where the quantity A is defined as 


A = 6n/e = 2(AV./ V.)/(J/6)(R/10)? (5.3) 


These figures indicate the following: 

(1) Increasing & from 0° to 90° generally increases 
fand f;; 

(2) Curves of f (or f,) vs. (& — &) are symmetrical 
about ( — &) = 180° only for the three cases & = 0°, 
90°, andz = 0°; 

(3) For 7 = 0, % is always an apo-apsis if either & = 
0° or i = 0°; ro is a peri-apsis for 7 = 90°; 

(4) For 7 = O, the number of apses per revolution 
can be two, three, or four depending on inclination 
and &p. 

Fig. 9 presents the relations between 7, €, 7, and & to 
assure an apse at (£* — &). The relationships between 
these variables are given by Eq. (4.24). This curve 
again illustrates the fact that an apse can occur at 
(— — &) = 180° only for & = 0° or 90°. 

Deviation from plane motion is illustrated in Fig. 10 
for the case & = 0°. The dimensionless quantity NV is 
defined by 

? o— (4/2) - 
Ns - - (5.4) 
(J/6)(R/ro)? sin 22 
where 6 is defined in Table | by Eq. (3.41). Positive 
deviations from the initial plane of motion occur for 
— from 0° to about 260° and negative deviations from 
260° to 360°, with a net negative deviation per revolu- 


tion. 
The injection speed required at a perigee altitude of 
300 nm above the equator (& = 0°) for various 


apogee altitudes, given by Eq. (4.27), is presented in 
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Fic. 10. Variations in the polar coordinate, 0, as a function of £ 
for several values of n (& = 0°). 
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Fic. 11. Speed required at a perigee altitude of 300 nm for a 
specified apogee (£ = 0) 


Fig. 11 for three cases: a polar orbit, an equatorial 
orbit, and the inverse-square force field. For a given 
’, — Y% (where rz is the radius at £ = 7m, a special case 
of r,) the injection speed is greatest for an equatorial 


orbit. 
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The Vertical Water-Exit and -Entry 
of Slender Symmetric Bodies’ 


JOHN P. MORAN* 
THERA, Incorporated 


Summary 


The hydrodynamics of a slender symmetric body in uniform 
vertical motion through the free surface of the water has been 
a linearized theory. Both thickness and lifting 
In the latter case, the analysis is applicable 


studied with 
problems are treated 
only to axisymmetric bodies; in the former, plane symmetric 
bodies are also treated. The analysis being based on successive 
approximations to the infinite-Froude-number situation, the 
results are accurate only when the Froude number is large. 
After the solutions are shown to consist of vertical distributions 
of singularities, applications of the theory to specific problems are 
illustrated. These include the determination of the time vari- 
ation of the forces and moments felt by a given body crossing 
Also 
exunined with the help of illustrative examples are the effects of 
the principal 


the surface, and a criterion for the avoidance of cavitation. 


the finiteness of the Froude number, which are 


contributions of the analvsis 


Symbols 
Ci upward thrust coefficient; see Eqs. (44) and (45) 
Cp = pressure coefficient; see Eq. (39) 
Cor see Eq. (41) 
f = doublet strength; see Eq. (22) 
WF integrals of f; see Eqs. (38) 
I upward thrust 
Fr Froude number U?/gl 
g = acceleration due to gravity 
l body length 
N(y) see Eq. (34) 
p = pressure 
Ps pressure on water surface 
q = pitch rate 
Q = source strength 
S(yv) = cross-sectional area = 7rX%(y) 
t = time 
U = body speed 
V(y) see Eqs. (30) and (36) 
x, ¥,@ = coordinates; see Figs. 1 and 2 
X(y) = body profile shape 
Ye = point about which g is measured 
Z = displacement of topmost point of body above un- 
disturbed position of surface 
a = angle of attack 
n* = location of single source 
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p density of water 
o = velocity potential 
T thickness-chord ratio 
w = frequency of pitching oscillations 
Subscripts 
6. 1, 2, ,n_ order in expansion of hydrodynamic effects in 


powers of Fr7! 


B due to buoyancy 

cr critical condition at which cavitation is in 
cipient 

W, A parts of infinite-Froude-number solution due 


singularities on body axis and on image 
body, respectively 
jee associated with thickness and lifting prob- 


lems, respectively 


Other Notation 


. in dimensional, space-fixed coordinates (Cor- 
responding unstarred quantities are in 
dimensionless body-fixed coordinates. ) 

7 derivative with respect to argument 

D{ j|/Dt convective derivative; see Eq. (5) 

O[ | order of | | 


Introduction 


pees OF THE HYDRODYNAMIC LOADS felt by a body 
striking a water surface date back to 1929,? their 
subsequent history having been reviewed by Szebe- 
hely.* More recently, the development of underwater- 
launched missiles has stimulated research into the re- 
lated water-exit problem.*~* 

In these previous studies, the effects of gravity on 
the fluid motion were ignored. That is, the Froude 
number Fr, the fundamental parameter of all problems 
involving an air—water interface, was assumed to be 
infinite. But in certain related problems which have 
been solved for large-but-finite Froude numbers, the 
effects of this parameter have been significant. For 
example, the slope of the lift curve for a flat-plate 
hydrofoil planing on the free surface was found to de- 
crease by 18 percent in going from the infinite-Froude- 
number case to one in which Fr = 20, and by a further 
9 percent when Fr = 10.*°' To fix ideas, note that 
if a body the size of Polaris were to leave the water at 
100 ft/sec, the Froude number based on its length 
would be about 10. 

Besides being restricted to infinite Froude numbers, 
most of the previous work is not applicable for arbitrary 
locations of the body with respect to the surface. the 
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a” Far from the body, ¢* must reduce to a constant. In 
Us addition, if the density of the air above the water sur 
face is neglected relative to that of the water below, ¢* 
0 must be such that the pressure is constant on the sur 


a 
x 






FREE 
SURFACE 
X(y) 


Vy 


Fic. 1. Geometry of thickness problem. O,O* origins of 
body-fixed and space-fixed coordinates, respectively. Y,x,y,Z 
dimensionless, normalized by / 


exceptions being the analyses of Goodman,* and Shiff- 
man and Spencer.'' That is, it is only useful either 
during the initial phase of the impact or while the body 
is completely submerged. 

It is the purpose of the present analysis to partially 
fill in these gaps in our knowledge of water-entry and 
-exit. Specifically, the problem of a slender symmetric 
body passing vertically through the free surface of the 
water will be treated for large, but finite, Froude num 


bers. 


(1) Formulation of the Problem 


(a) Scope of the Study 

The physical situations to be studied are depicted 
in Figs. 1 and 2. The former figure shows a plane, 
symmetric airfoil or a body of revolution in vertical 
constant-speed passage through a horizontal water 
surface. In the latter, a body of revolution is shown 
executing some lateral and pitching motions as it crosses 
the surface. In parlance familiar to aerodynamicists, 
both the thickness and /ifting problems of water-exit 
and -entry will be treated for axisymmetric bodies, but 
only the former problem for the symmetric airfoil. 

To permit a unified discussion of these problems, we 
define a space-fixed coordinate system (x*, y*, 0*) as 
follows: Let y* be measured vertically upward from 
the undisturbed position of the free surface. In the 
plane flow case, (x*, y*) is a Cartesian coordinate sys- 
tem and 6* is meaningless. In the three-dimensional 
situation, x* is measured radially and 6* is a meridional 
coordinate measured from the plane of the motion. 


(0) Equations, Assumptions, and Boundary Conditions 

We assume that the flow is incompressible and irrota- 
tional. Thus we may work in terms of a velocity po- 
tential @* which, from continuity, satisfies Laplace's 
equation 


Ve" =. (1) 


face. Assuming the body is slender enough so that 
the fluid motion represents a small perturbation of a 


state of rest, we may write this condition as follows:'? 


076* /Ot*? + 2 06*/Ov* = 0 (2) 


Consistent with the small-perturbation assumption, 
Eq. (2) is applied at y* 0 rather than on the free 
surface itself. 

It is somewhat more convenient to work in the di 
mensionless body-fixed coordinates (x, y, 0; Z) which 
as shown in Figs. 1 and 2, are related to (x*, y*, 6*: 
t*) by 

f= g7 I, y (Z — y*) /i, 6 = 6%, (3 
Z= +Ut*/I, @ = ¢*/Ul 


Here / is the body length, L’ its speed, and Z the dis 


- 
a 
i 
. | 
an 
x ’ 2 
- “es 
~ 
4 X(y) 
‘y 


Fic. 2. Geometry of lifting problem. All coordinates di- 


mensionless, normalized by /. Body moving in x,y plane, undis 
turbed position of free surface lies in x). plane 


placement of its topmost point above the unaisturbed 
position of the free surface, the upper sign being ap- 
plicable to the exit problem and the lower to the entry 
situation. 

In these coordinates, Eq. (2), the free-surface bound 


ary condition, becomes 
D*¢/Dt? — Fr- d/dy = 0 [fy = Z] (4 
) Dt is a convective derivative 


)/O0Z (9 


where D( 

Dt )/ Pt oC ) oy + A 

and Fr is the Froude number (or, according to some 
definitions, its square), 

Fr = U?/gl (6) 

We define Y(y) as the dimensionless ordinate of the 

airfoil or radius of the body of revolution, as the case 

may be. The linearized condition of no flow through 
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the body may be written in a unified form valid for both 
cases as 
| cos 6 

ix = X(y)] (7) 


0p Ox = X'(v) — [a + HY — V2) 


Here ais the angle of attack of the body and g its pitch 
rate about y = y,, both quantities being zero in the 
plane flow problem by definition of the scope of the 
study. For convenience, we have supposed the lateral 
and pitching motions of the body of revolution to be 
coplanar, but this will not alter our conclusions. 


(c) Equivalence of Exit and Entry Problems 


We assume that the flow is noncavitating, and 
thereby avoid the complication of specifying an addi- 
tional boundary condition (of constant pressure) on the 
free surface of the cavity. Further, coupled with the 
assumption that the fluid is inviscid, this restriction 
enables us to equate the exit and entry problems. 
That is, if the vertical component of the body’s velocity 
is reversed, but all other conditions held constant (in- 
cluding the orientation of the body with respect to the 
surface and the lateral component of the body's mo- 
tion), the net force felt by the body is unchanged. '* 
Thus, although this paper is applicable to problems of 
both water-entry and -exit, we are able to simplify 
our derivations by referring them solely to the exit 
problem. We have anticipated this in the writing ot 
Eq. (7), which is correct only for an exit situation. 


(d) Separation of Thickness and Lifting Problems 

In view of the linearity of Laplace’s equation and the 
boundary conditions, we may split the potential @ 
into two parts, @7 and @¢, (say), such that both parts 
satisfy Laplace’s equation and the free-surface bound- 


ary condition, Eqs. (1) and (4). On the body, how- 


ever, 

07 /Ox = X’(y) [x X(y) | (8) 
while 
Og, ‘Ox —la+q(y—y.)|cosé [x = X(y)] (9) 


The potential ¢7, which from Eq. (8) is clearly sym- 
metric about the body axis x = 0, may be determined 
in an unbounded-fluid problem by distributing sources 
along that line. Such a determination constitutes the 
thickness problem of slender-body theory, while the 
antisymmetric ¢, is associated with the /ifting problem, 
which is conventionally solved by a distribution of vor- 
tices (in two dimensions) or doublets (in three dimen- 


sions). 


(e) Delineation of Froude Number Effects 

From its definition, Eq. (6), the Froude number is 
seen to be the ratio of the free-stream dynamic pres- 
sure to the gravity-induced pressure differential from 
one end of the body to the other. Thus, it may be 
regarded as a measure of the importance of hydro- 
dynamic effects relative to hydrostatic effects. The 
latter effects, although of prime importance in problems 


involving surface vessels, exert a relatively minor 
influence in problems of water-exit and -entry. It 
thus seems sufficient to seek an approximate solution 
of the problem which is valid for large, but finite, Froude 
numbers. 

To this end, we assume that the potential is expand- 
able in a power series in the reciprocal of the Froude 
number, the leading term of which is the infinite- 
Froude-number solution. 


r+) pv o, Fr ; (10) 

n=0 
Here each of the ¢,'s must be a solution of Laplace’s 
equation, due to its homogeneity and linearity. To 
obtain free-surface boundary conditions on the @¢,’s, 
we substitute the assumed form of the solution into 
Eq. (4) and equate terms of like order in Fr~!. We 


thus require 


D*o Dt? () [y Z| (11) 
D*¢, Dt* Op) OV |) Z| (12) 
D*o, Dt? O¢,,1 OV Lv z.n> 1 (13) 


(2) The Infinite- Froude-Number Solution 


Eq. (11), the free-surface boundary condition on the 
infinite-Froude-number solution, is readily reduced to 
the requirement that ¢» (or @*)) be constant along the 
y = Z (or y* = 0) plane. The solution of such a prob- 
lem is well known to be expressible in terms of anti 


symmetric images. 


(a) Single Singularity Rising Toward Surface 

For example, suppose we have a source of strength 
Q in motion below the surface, as shown in Fig. 3. The 
free-surface condition then may be satisfied by posi- 
tioning a sink above the surface so that the y* 0 
plane is the perpendicular bisector of the line joining 
the two singularities. For, as indicated in Fig. 3, the 
resultant velocity induced by the source-sink combina- 


tion is purely vertical at y* = 0, so that 0¢*)/Ox* = 0 
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Fic. 3. Infinite-Froude-number solution to problem of single 
source in motion below surface 
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and ¢*) = constant on that plane. Thus ¢*) may be 
split up into two parts, one being due to the singularity 
below the surface and the other to its image. Desig- 
nating these parts by $*w and $*4, respectively, we 
have for the problem depicted in Fig. 3, 


o*w = (Q/4m) In [x*2 + (y* — n*)?] (14) 
o*, = —(Q/4zm) In [x*2? + (y* + 9*)?] = (15) 


if the source is two-dimensional and located at (0, *). 

Similar results may be worked out for the situation 
where the source is three-dimensional, to which case 
the above discussion is equally applicable. Further, 
if the source is replaced by a doublet, the free-surface 
condition is readily shown to be satisfied by replacing 
the sink by an oppositely directed doublet, while two 
vortices of the same sign also satisfy Eq. (11) on the y = 


Z plane. 


(b) Solid Body Rising Toward Surface: 
Thickness Problem 

Now suppose that instead of a single singularity, we 
have a symmetric body in vertical motion toward the 
surface. Insofar as the thickness problem is concerned, 
we anticipate a distribution of sources along the body 
axis. For each member of the distribution, we must 
position a sink at its image point in order to satisfy 
Eq. (11), so that the net result is a distribution of 
sources and sinks along a vertical line. But if the 
body is sufficiently slender, the boundary condition 
of no flow through the body at a given point thereon 
influences the source strength only at that point. 
Thus, only the sources below the surface—i.e., on the 
bod y—enter into the satisfaction of this boundary condi- 
tion and so are of a strength which is the same as if the 
fluid were unbounded. Physically, the flow due to 
the singularities above the surface is directed mainly 
parallel to the actual body and, thus, within the linear- 
ized approximation, has a negligible velocity com- 
ponent normal to the body surface. 

To express this in mathematical form, we distribute 
sources whose potential is given by Eq. (14) from n* = 
Z* to n* = Z* — |, apply the transformation to di- 
mensionless body-fixed coordinates, and note that in 
an unbounded fluid the source strength is 2X’(y), 
thereby obtaining 


ya 
orw = on f X'(n) In [x? + (vy — n)*J dn (16) 


Similarly treating Eq. (15), we find the zeroth-order 
correction to ¢7y due to the presence of the free surface: 


ew 
or, = — | X'(m) In [x? + (vy + 9 — 2Z)?] dn 
2r J0 
(17) 


The corresponding formulas for the axisymmetric case 
are as follows: 


1 (' 
ory = — i} S’(n) [x? + (vy — 9)?]—1/2 dy (18) 
4nr JO 
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1 
or, = { S’(n)[x? + (vy + 9 — 2Z)?]-/2 dy (19) 
4r J0 
where S(y) is the dimensionless cross-sectional area of 


the body, rX*(y). 


(c) Solid Body Rising Toward Surface: Lifting Problem 

In the two-dimensional lifting problem, we again 
anticipate a distribution of singularities along the chord 
of the airfoil. But here the situation is readily seen to be 
more complex, since the image vortices induce a flow 
which has a nonnegligible component normal to the 
body surface. Thus we must also anticipate a ‘“‘wake’”’ 
behind the airfoil, the vortices composing which would 
be used to help satisfy the body boundary condition, 
and see further that there will also be a ‘“‘wake’’ behind 
the image airfoil. 

The lifting problem for a body of revolution is much 
simpler, however. For, just as in the source distribu- 
tion problem, slender-body theory tells us that the 
boundary condition at a given point on the body in 
fluences the doublet strength only at that point. Thus 
the image doublets may be ignored in satisfying Eq. 
(9), and the infinite-Froude-number solution may im 
mediately be written down as follows: 


I a 
ory = x cos 6 [ f(n) [x2 + (y — n)?]—3/2 dn 
tor J0 


(20) 


Pe 
or, = ’ x cos 6 { f(n) (x? + (y + 9 — 2Z)?]8/? dy 
4 Jo 


(21) 
The doublet strength f(y) is the same as though the fluid 
were unbounded, viz., 





f(y) = 2[a + qv — »-) |S(y) (22) 


. 


(d) A Conflict in Boundary Conditions 

Before moving on to the effects of the finiteness of 
the Froude number, it is convenient to use the simpler 
infinite-Froude-number solution to point out a failing 
of the linearized approach. For in this case, the con- 
stant-potential free-surface boundary condition is 
readily seen to conflict with the body boundary condi- 
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Fic. 4. Pressure distributions on biconvex airfoil in vertical mo- 
tion below water surface. 7 = 0.05, Fr = 10. 
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tion, Eq. (7), while the body is broaching the surface 
(Z > ©). 

Before we attempt to support the acceptance of this 
situation, a few words are in order regarding the form 
of the solution for Z > 0, to which case Eqs. (16) 
through (21) are inapplicable. The portion of the 
body above the surface, being in a relative vacuum, 
clearly cannot affect the flow below. Thus, of the 
singularities distributed along the body axis, we must 
count only those below the surface. The infinite- 
Froude-number solution then consists of the submerged 
portion of the body singularity distribution joined at 
the y = Z plane to its antisymmetric image. 

Now, as noted previously, in each of the problems 
treated in this paper, the body singularity strength is 
determined by satisfaction of only the local body 


boundary condition. But at y = Z, the body singu- 
larity is cancelled by its image, so that Eq. (7) cannot 
be satisfied at this point. At y = Z + e, however, 


for ¢ arbitrarily small, the strength of the body-bound 
singularities is unaffected by their images, so that there 
the body boundary condition can be satisfied. Thus 
the conflict in the boundary conditions is confined 
to the point y = Z. 

These statements are accurate within the linearized 
approximation—i.e., to the order of s—where 7 is the 
body thickness ratio. Clearly, the true extent of the 
region in which our solution is in error is also of order r. 

This fact has been utilized in reference 14 to estimate 
the errors induced by this conflict in boundary condi- 
tions. The pressure distribution predicted by the 
slender-body theory was modified in the region between 
y = Zand y = Z + 7 to agree with the constant- 
pressure free-surface boundary condition. The time 
variations of the lateral force and pitching moment 
were then calculated from the two pressure formulas. 
Although the modified pressure distribution differed 
quite radically from the theoretical one, the net effect 
in the force~and—moment-vs.—time curves was a small 
shift in the time scale. Thus it is felt that the violation 
of the body boundary condition at y = Z will not greatly 
impair the accuracy of our results. 


(3) Extension of Solution to Large-but- Finite 
Froude Numbers 


Now from Fig. 3, or from Eqs. (16)—(21), it is seen 
that the vertical components or the velocities induced 
by the singularity of the body below the surface and 
by its image are equal at the surface. Thus 

Odw - Og. (23) 
oy y=Z Ov y=Z 
and Eq. (12) may be rewritten as 
D? fe) ; 
$1 _ 5 Oba (24) 


Dt? y=Z + oy y=Z 


The techniques used to determine ¢; and higher order 
contributions to Eq. (10) will be illustrated in detail 
for the two-dimensional thickness problem. After- 
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Fic. 5. Pressure distributions on biconvex airfoil in vertical mo 
tion through water surface. +r = 0.05, Fr = 10 


ward, the results for the problems involving bodies of 
revolution will simply be written down. 


(a) Plane Thickness Problem 


We first substitute the assumed form of the solution, 
Eq. (10), into the body boundary condition, Eq. (8). 
The latter relation is inhomogeneous in @ and does not 
involve Fr, and so 


O¢r,/Ox = 0) [x = X(y), n>O] (25) 


Now due to the nature of the problem, the flow asso- 
ciated with the ¢7,’s must be symmetric about the x 
axis. The singularities inducing this flow must, there- 
fore, be sourcelike—i.e., if not actually sources, they 
must be doublets (or higher multipoles) directed along 
the y axis. But a distribution of such singularities 
along a line creates a discontinuity in the velocity nor- 
mal to that line, and therefore, from Eq. (25), the 
orn's (n > O) cannot be singular on the body axis. It 
is then evident that just as in the infinite-Froude- 
number case, the strength of the sources distributed 
along the airfoil axis is the same as if the fluid were 
unbounded. 

The ¢7,’s may be regarded as the real parts of com- 
plex potentials, so that Eqs. (13) and (24) are relations 
between analytic functions of a complex variable which 
are valid along a plane. Therefore, these relations may 
be continued off the y = Z plane, becoming differential 
equations rather than boundary conditions. Note that 
if Eq. (23) were used to eliminate $7, from Eq. (24) 
in favor of dry, the resultant boundary condition could 
not be so continued. For @7, is, by definition, singular 
on the body axis, on which such quantities of interest 
as the pressure are computed in the linearized theory. 

We now substitute Eq. (17) into Eq. (24) and obtain, 
after an integration by parts, 

Dor, 2 [ P x? — (y+ — 2Z)? 7 
Dt =x Jo X(n) [x* + (y + 9 — 22Z)?}?  .o 


This relation bears a marked resemblance to D*gr, ‘D??, 
which from Eqs. (5) and (17) may be written, 


D°or4 l . os x? — 7+ eo 2Z)? 
a X'(n) -, —— ay 

Dt? rJ0 [x2 + (y + 9 — 2Z)?)? 
(27) 
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Noting that the only differences between Eqs. (26) and 
(27) are a factor of (—2) and an integration of the 
singularity strength, we similarly modify Eq. (17) and 
obtain for @7,. 


ze 
or, = | X(n) In [x2 + (vy + 9 — 2Z)?]dn (28) 
T « 


0 


We then substitute this result into the differential 
equation for ¢7,, which is obtained from Eq. (13), and 


again integrate by parts. 


Don _ 2) §+1~ 

De x [x? + (y + 1 — 2Z)?] 
9 £3 , eee — 97)2 
= | Vin) x ye a = ) dn (29) 
rJ0 [v2 + (y + n — 22Z)?/? 


Here, 


ro) = | X(¢) dé (30) 
J 0 


The integral part of Eq. (29) may be satisfied by in- 
cluding in 7, a term, 7,” (say), which is of the same 
form as Eqs. (17) and (28); viz., 


ve 
¢r(! = V(m) In[x? + (y + n — 2Z)?] dy 
res/0O 


(31) 


The nonintegral term of Eq. (29) may be integrated to 
yield for the remaining part of @7,, ¢7,'°’ (say) : 


227-1 
or? = — Vil) | In (x? + (y — ¢)?] de (32) 
T Jz 
Here we have introduced a change of variable so that 
gr.” is recognizable as the potential of a source dis- 
tribution of constant strength running from the image 
point of the tail of the airfoil (vy = 2Z — 1) to the free 
surface (y = Z). 

We could have let this source distribution run from 
y = 2Z — 1 to y = —o, which would also satisfy 
Eq. (29). However, such a distribution would grow 
in extent as the body left the water, whereas the poten- 
tial given by Eq. (32) disappears as Z — 1, which is 
physically more plausible. * 

The solution given above is valid only while the air- 
foil is completely submerged. For as previously noted, 
the members of the body source distribution are ‘‘ex- 
tinguished’’ as they pass through the surface and simi- 
larly for those singularities distributed over the image 
body. Egs. (16), (17), (28), and (31) are therefore 
modified in that the integrals now extend from n = Z 
(the free surface) ton = 1. Substituting the resultant 
expressions back into Eq. (4), we find that it is also 
necessary to add a source of time-dependent strength 


* Also, in this two-dimensional problem, it is found that the 
alternative solution violates the boundary condition that @7 
reduce to a constant far from the body. In the axisymmetric 
case, however, both versions of the solution satisfy all the 
boundary conditions of the problem, and the nonmathematical 
reasoning given above must be used to choose between them. 
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at the intersection of the airfoil trajectory and the un 
disturbed position of the free surface. We thus obtain 
the final result as 


St 
7 X'(n) In [x? + (y — n)?] dy — 
2a J1(Z)Z 


Bw 
z | [.X'(n) — 2Fr-! X(n) + 2Fr-°V(n)] X 
Te 


1(Z)Z 


In [x? + (vy + » — 2Z)?| dn — 


227-1 
Fr-? V(1) in [x*? + (y — £)71 df + 
T JZ 


1(Z) — | Fr-! V(Z) — Fr-? N(Z)} X 


In [v? + (y — Z)?] + O[Fr-] (33) 


where Ny) = | V(¢) dé (34 
/J0 


Eq. (33) has been made valid for the two cases of com 


plete submergence and broach by introduction of the 


unit step function, 1(Z). 


(b) Problems Involving Bodies of Revolution 

Now turning to the axisymmetric case, we assume 
that the free-surface boundary conditions on the 
various ¢,'s, Eqs. (13) and (24), may be continued into 
the flow field below the surface, just as in the two 
dimensional case. Here the continuation is not 
analytic; our justification is a posteriori, in that the 
solution obtained satisfies Laplace’s equation and 
all the boundary conditions of the problem. 

We first seek the axisymmetric analog of Eq. (33), 
the solution to the thickness problem. In the two- 
dimensional case, we implicitly assumed X(1) = 0—1.e., 
that the airfoil was closed. From a “‘practical’’ point 
of view, such a stipulation is not overly restrictive. 
However, one often encounters bodies of revolution for 
which X(1) # 0. An admittedly approximate pro- 
cedure for handling such situations is to tack onto the 
tail of the body an infinite circular cylinder of the base 
dimensions. Such a step is clearly most nearly valid 
when the slope of the body is zero at the tail. 

It this approximation is made, by following through 
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Fic. 6. Maximum speed of biconvex airfoil crossing water 
surface vertically without cavitation. pp, = 14.7, 0.2 psi, 
respectively. 
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Fic. 7. Upward thrust on biconvex airfoil crossing water surface 
vertically. + 0.05 


the same procedure detailed above for the two-dimen- 
sional case, we may find an expression for the potential 
about a slender body of revolution which is valid to 
the order of Fr~' if the body is open-ended [X(1) # 0} 


and to the order of Fr~? if it is closed [X(1) = 0}: 
z 

+ = — | S’(n)[x? + (vy — 9)?]7/?2 dn + 
4a J1(Z)Z 


By 
| [S’(n) — 2Fr-! S(n) + 2Fr-2V(n)] X 
we i(Z)Z 


l 
[x2 + (y + 9 —2Z)?]—"/2 dn — 2 [Fr-? S(1) — 
T 


22Z—1 
Fr-2 V( 1)] [x? + (y—- ¢)?] l 2 de a 
JZ 


| 
1(Z) [Fr-! V(Z) — Fr-? N(Z)] X 
or 


[x? + (y — Z)2]-¥2 


Here N(y) is as given by Eq. (34), but I’(y) is rede- 


Vy) = [ S(¢) dt (36) 
/7 0 


The solution of the corresponding lifting problem 
may be found by proceeding along similar lines. Here 
only a first-order correction to Eqs. (20) and (21) was 


fined as 


found, the result being 


l 


+7 


l 
x cos 6 / f(n) (x? + 
J1(Z)Z 
*] 


l 
(y — )?|~3/2 dyn + x cos 6 | [f(n) — 
i dr J 1(Z)Z 
2Fr-! f(y) |[x? + (y + n — 2Z)?]-3!? dn — 


227-1 
x cos 6Fr-! f(1) } [x*§ + Gy — fT et — 
e Z 


Le 
‘ ] 

1(Z) = x cos 6 X 
2r 


Fr-f(Z) [x2 + (y — Z)?]-3/2 + O[Fr-?] (37) 
“Fe ' y . 
where fo = | f(©dg, fy) = [ f(g) dé (38) 
0 /J0 


and f(y) is given by Eq. (22). 
In this result, it is implied that f(y) is time-inde- 


-EXIT 


AND -ENTEY SO9 


i.e., that the angle of attack a and the pitch 
If, instead, these quantities are 


pendent 
rate g are constant. 
functions of time, it is found that the infinite-Froude- 
number solution is unchanged, but the determination 
of $,, is considerably less straightforward. This has 
been carried out for the case where a and g are slowly- 
varying functions of time,'* as would be the case if the 
body were executing oscillations in its angle of attack 
with small reduced frequency w/ LU. 


(c) Summary of Solutions for Velocity Potential 


It is of interest to note that the three solutions, 
Eqs. (33), (85), and (37), each consist of the following 
elements: 

(1) a distribution of singularities (sources in the 
thickness problems and doublets in the lifting problem) 
along the axis of the submerged portion of the body 
the singularity strength being the same as if the fluid 
were unbounded; 

(2) a distribution of the same singularities along the 
axis of the image body, the strength of the distribution 
being a function of Froude number and body shape; 

(3) a constant-strength singularity distribution from 
the tail of the image body to the free surface*; and 

(4) while the body is broaching the surface, a singu- 
larity at the intersection of the body trajectory and 
the free surface, the strength being a function of both 
time and Froude number. 


(4) Applications of the Theory 
(a) Pressure Distributions 
We define a pressure coefficient, 
C, = (p — ps)/(1/2pU”) 39) 


where /p, is the constant pressure along the water sur- 
face. With the knowledge of the potential, the C, dis- 
tribution over a given body may be calculated by em 
Bernoulli’s 


ploying suitably linearized versions of 


equation. In the plane, symmetric case, we have 


Cp = 2Fr-' (vy — Z) — 2D¢ Di (40) 


with the first term representing the hydrostatic pres- 
sure on the body. For the axisymmetric thickness 
problem, D@/ Dt is of the same order in 7 near the body 
as is (0@ Ox)*. Thus in that case, it is necessary to in 
clude in Eq. (40) an additional term, 
>, = — (0¢/Ox)? (41) 
Similar quadratic terms must be added to the right- 
hand side of Eq. (40) in the lifting problem for an 
axisymmetric body. 

Example: In Figs. 4 and 5, we have plotted some 
typical pressure distributions computed from Eqs. 


* This is highly reminiscent of the solutions obtained by 
Havelock and Tan*® for the flow about various submerged 
bodies and singularities moving parallel to the surface. In 
those studies, it was found necessary to introduce a horizontal 
stream of singularities trailing behind the image body 








810 JOURNAL OF THE 


(33) and (40) for the case of a biconvex airfoil formed 
by two parabolic arcs, the profile shape of which is 
given by 

X(y) = 2ry(1 — y) (42) 


In the submerged case (Z < 0), the pressure distribution 
is not noticeably affected by the presence of the free 
surface (except for a shift in level due to hydrostatic 
effects) until the airfoil is very close thereto—for the 
case shown, until it is within about a tenth of the air- 
foil chord. The singular behavior of C, at y = Z in 
the broach case (Z > () reflects the conflict in the 
boundary conditions at that point noted previously. 
For y < Z, C, is zero, by assumption. 


(b) Inception of Cavitation 

In this analysis, it has been assumed that no cavita- 
tion is present. It is therefore important to have a 
criterion for the inception of cavitation, so as to estab- 
lish limits for the applicability of the theory. Also, 
in view of the well-known drag increases which ac- 
company cavitation," such a criterion could be used to 
define operating limits for the vehicle under considera- 
tion. 

For our present purposes—to illustrate possible 
applications of the theory—it is sufficient to ascribe 
the appearance of this phenomenon to the lowering 
of the fluid pressure below the vapor pressure. The 
maximum speed of a given body above which it will 
cavitate, U,,, may then be determined as follows: 
From data such as that shown in Figs. 4 and 5, obtain 
the minimum pressure coefficient felt by the body as 
it crosses the surface. This will be a function of 7 and 
Fr. Set the pressure in the definition of C, equal to 
the vapor pressure, p,. From Eq. (39), we then ob- 
tain 

U.2(1, Fr) = 2 (bo — ps)/pComin(7, Fr) (48) 
The corresponding length of the body may then be 


determined from the definition of Fr, Eq. (6). We 
thus obtain, though not directly, U,,(7, /). 























Example: In Fig. 6, the critical speed of the bicon- 
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vex airfoil has been plotted as a function of its length 
for two thickness ratios, p, and p, being taken as ().2 
and 14.7 psi, respectively. Since the calculations are 
based on a theory which assumes the Froude number 
to be large, we have also included in the figure curves 
of constant Froude number. 

It was found that, to a very good approximation, 
U., ~ 7r'*. From Eq. (43), this indicates that 
Comin 
dynamic contribution to the pressure would vary 
with 7, it is remarkable that the inclusion of the hydro- 
static pressure does not ostensibly alter this dependence 
oi Cy. ,,, OF. 

It should also be noted that for this particular body, 
the critical condition occurred when the airfoil nose 
In general, however, 


~ r. Although we would expect that the hydro- 


just touched the surface (Z = 0). 
the value of Z for which C, is a minimum is a function 


of Fr and r. 


(c) Upward Thrust 
We define an axial thrust coefficient, Cy, as the ratio 
of the upward thrust F felt by the body to the product 
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Fic.9. Force functions for biconvex airfoil crossing water surface 
vertically. See Eq. (46). 


of its maximum frontal area and the free-stream dy- 
Thus in the plane case, 
F 


C= = (44) 
°* (pU?2/2) al 


namic pressure. 


while for the body of revolution, 


i al , 

* (pU2/2)(w72l2/4) _ 

Ordinarily, Cy may be determined by a straight- 
forward integration of the C, distribution. However, 
for certain bodies of revolution—specifically, for any 
body with a nose or tail as blunt as or blunter than 
the linearized approach breaks 
The origin of this 


that of an ellipsoid 
down, predicting an infinite thrust. 
failure is Cre which becomes too strongly infinite 
near the blunt end(s) of the body to integrate. 

In Fig. 7, the time variation of Cp is 
Far below the sur- 


Examples: 
shown for the biconvex airfoil. 
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Fic. 10. Force functions for body formed by revolving a para- 
bolic are crossing water surface vertically. Solid lines, 7 = 0.10; 
dotted, r = 0.05. See Eq. (46 


face, the airfoil feels only a buoyant force, which in 
dimensionless form is proportional to Fr~!. The in- 
crease in Cy as the airfoil approaches the surface may 
be ascribed to a gradual transfer of kinetic energy from 
the fluid to the body. Alternatively, observe that the 
fluid in front of the body becomes less ‘“‘stiff’’ as Z > 0, 
since it is backed up by a relative vacuum. 

The discontinuity in the slope of Cy at Z = 0 may 
be rationalized by a consideration of the effects of the 
singularity-destruction process which takes place for 
Z> 0. Specifically, consider the contribution of the 
sources on the body. While the body is completely 
submerged, from symmetry these sources do not con- 
tribute to the axial thrust. On entering the broach 
phase of the motion, however, a source is immediately 
destroyed, so that the remaining sources on the body 
suddenly contribute positively to Cr. <A similar proc- 
ess is simultaneously taking place with the image 
sinks of the infinite-Froude-number solution, the effect 
of which is in a direction opposite to that of the sources 
on the body, and, from Fig. 7, is apparently stronger. 

For Z > 0, the shape of the thrust curve is not un- 
like that of a roller coaster at the larger Froude num- 
bers, which may also be rationalized by a consideration 
of the singularity-destruction process. The smoothing 
out of this shape at the smaller Froude numbers is 
largely due to the increased prominence of the buoyant 
force, which decreases monotonically as the body 
leaves the water and, as noted previously, is propor- 
tional to Fr~!. 

The time variation of the upward thrust felt by the 
body formed by revolving a parabolic arc, whose pro- 
file shape is the same as that given by Eq. (42), is 
shown in Fig. 8. Except at Z = 0, the shapes of the 
thrust curves in Figs. 7 and 8 are quite similar. That 
Cy’(Z) is now continuous at Z = 0 is due to the fact 
that the cross-sectional area distribution of this body 
is cusped at the nose, so that the source strength is 
zero at that point. Hence the broach phase of the 
motion is entered quite smoothly, without the immedi- 
ate destruction of a singularity as in the previous 


case. 


Note that in view of the equivalence principle quoted 
earlier, Figs. 7 and 8 are equally applicable to the 
water-entry problem. If the bodies studied had 
lacked fore-and-aft symmetry, however, we would have 
required two sets of curves for each orientation of the 
body—1.e., nose up and tail up. 


(d) Lateral Force and Pitching Moment 


Let us now turn to the application of Eq. (37) to the 
lifting problem. It will be recalled that our solution 
to this problem consisted entirely of a distribution of 
doublets along a vertical line, with those doublets on 
the body being of the same strength as if the fluid were 
unbounded. But from slender-body theory, we know 
that the difference in pressure across a lifting body at a 
given point thereon is determined by only the local 
singularity strength, provided that strength is time- 
independent. Since this result must hold as well in 
the present case, we conclude that the lateral force 
felt on a given section of the body is independent of 
both its distance from the free surface and the Froude 
number. The net lateral force and the pitching mo- 
ment will thus be the same as if the fluid were un- 
bounded while the body is submerged, but will decrease 
to zero as the body leaves the water. This same re- 
sult was obtained by Goodman,* who analyzed the sta- 
bility derivatives of a slender body vertically crossing 
the surface at infinite Froude number. It is remark- 
able that this conclusion is unchanged by the finiteness 
of the Froude number. 


(e) Effects of Finiteness of Froude Number 


In order to pin down the effects of the finiteness of 
the Froude number, the expression for Cy was broken 


up to read, 
: — T ae 
Cr = Crs Fr + | |  B Cy, Fro" (46) 
LT Jj »=0 Z 


where the upper term in the square brackets applies 
to the two-dimensional situation and the lower to the 
axisymmetric case. Here C,y,/r~' is the buoyant 
force. In Figs. 9 and 10 the various force functions 
Cr, are plotted against time for the biconvex airfoil 
and for the revolved parabolic arc, respectively. In 
the former case, the force functions are completely 
independent of both 7 and Fr, while in the latter, there 
is a weak dependence on 7 due to the appearance 
of terms proportional to In 7 in the C,,’s. 

For the airfoil, all the C,,,’s are generally of the same 
order of magnitude, while for the body of revolution, 
Cy, is dominant. Since 7 and Fr~! must be small for 
the theory to be applicable, it is thus apparent that a 
‘modified infinite-Froude-number approximation” for 
Cr, 

. _— , i ~ 
Cr ~ Cr;Fr '+[%.] ce (47) 
{7 
would be almost as accurate as Eq. (46), particularly 
since C,y,’s are available only for n <2. 
Although numerical comparison of the two formulas 
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tends to substantiate this supposition, it should be 
pointed out that the effects of the finiteness of the 
Froude number can be important in establishing certain 
details of the pressure distribution. For both the 
biconvex airfoil and the revolved parabolic arc, the 
“modified —infinite-Froude-number — approximation” 
yielded errors of about 8 percent in the level of the 
minimum pressure coefficient at Fr = 10, at least as 
compared with the values calculated to the order of 
Fr~*. The importance of such a quantity has been 
noted previously in connection with the cavitation- 
inception criterion. 

The possibility that Eq. (46) may be readily trunc- 
ated is quite encouraging, since it is indicative of rapid 
convergence of the large-Froude-number expansion. 
Since only the first few terms of this expansion were 
found, convergence could not be studied mathemati- 
cally. 


Conclusions 


The solutions obtained depended on the following 
three steps: 

(1) The body under study was assumed to be slender 
justify a complete linearization of the 
localized 


enough to 
boundary conditions. 
conflict in these boundary conditions, but the effects 
of this violation are felt to be small. 

(2) The Froude number was supposed to be large 
enough to justify the expansion of the solution in powers 
Only the first few terms of the ex- 


This resulted in a 


of its reciprocal. 
pansion were found. 

(3) After the influence of the body-bound singulari- 
ties was removed from the velocity potential, the free- 
surface boundary condition was continued into the 
entire region below the surface. This continuation 
was analytic in the plane problem, but constituted an 
important assumption in the cases involving a body of 
revolution, whose justification was a posteriori. 

Using a_ singularity-superposition approach, we 
found that the solution consisted entirely of a dis- 
tribution of singularities along the axis of motion. 
Those beneath the surface (on the body) were of the 
same strength as if the fluid were unbounded. Those 
above described, in part, a distorted image body, 
whose shape was a function of Froude number. Also 
involved were a trail of singularities from the free 
surface to the tail of the image body, and a discrete 
singularity at the intersection of the surface and the 
path of the body’s motion. 

The principal contributions of 
effects of the finiteness of the Froude number 
found to be quite small in all cases examined, at least 
insofar as the net forces felt by the body were con- 
Indeed, within the linearized approximation, 


the 
were 


the analysis 


cerned. 
+ 
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these effects may be neglected entirely in calculating 
the lateral forces and pitching moment. However, 
the Froude number did exert a fairly strong influence 
on certain details of the pressure distribution, such 
as the minimum pressure coefficient. This quantity 
was important in calculating an approximate criterion 
for the avoidance of cavitation. 
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On a Tvpe of Clamped Edge Condition for 
Plate Problems’ 


F. ESSENBURG* anp R. C. KOELLER** 
Illinots Institute of Technology 


. 
Summary 

A method is given for the treatment of plates, portions of 

which are constrained by two smooth, rigid, plane, clamping sur 

faces. By employing the E. Reissner plate theory, expressions 

for the stresses and displacements, valid for the region constrained 


within the clamp, are obtained and appropriate continuity condi 


tions are satisfied at the edge of that region. The numerical re 


sults of a specific axisymmetric example illustrate the sub 


stantial improvement of the present treatment, 


(1) Introduction 


_ PRESENT paper is concerned with the bending of 
uniform, isotropic, elastic plates, portions of which 
are clamped between two smooth rigid plane surfaces. 
In particular, the analysis of the present paper is con- 
cerned with axisymmetric circular plates, having an- 
nular portions constrained by smooth, rigid, clamping 
members. Consistent with the assumption that the 
transverse normal component of the displacement of 
the plate is independent of the thickness coordinate, to- 
gether with the usual assumptions associated with 
smooth rigid surfaces, it will be assumed that those por- 
tions of the plate contained within the smooth clamp 
will satisfy the following conditions: (1) The deflection 
is constant; and (2) the top and bottom surfaces of the 
plate are free from shearing tractions. Such assumed 
conditions serve as an idealization for the clamped edge 
conditions for the particular case where the interior sur- 
faces of the clamping members are relatively smooth 
and where the rigidity of the plate is small in compari- 
son to the rigidity of the clamping member. While the 
usual clamped or built-in edge conditions! are physi- 
cally consistent (as well as experimentally verified) for 
the case where the surface of the clamp is sufficiently 
rough to prevent slippage, it should be noted that al- 
though they imply the application of stress couples and 
resultants at the edge of the clamp, no information is 
given regarding the displacement and stress fields in 
that portion of the plate contained within the clamp. 
For the case of a smooth clamp, the method of the 
present paper is to obtain expressions for the stresses 
and displacements, valid for the plate region contained 
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within the clamp, and to satisfy appropriate continuity 
requirements on the stresses and displacements at the 
edge of that region. The Poisson-Kirchhoff plate 
theory cannot be utilized in such an approach, for it will 
be recalled that, as a consequence of the restrictive as 
sumptions of the classical theory, the stresses and dis- 
placements are all related to the derivatives of the de- 
flection. Thus for the region of the plate contained 
within the smooth clamp, where the deflection is con- 
stant, the classical theory predicts that the stresses 
and displacements will all vanish. Continuity of such 
a (vanishing) displacement field at the edge of the 
clamped portion (where it must also be concluded that 
the stress vector field is discontinuous) yields a result 
identical to that obtained by the application of the usual 
clamped edge conditions, and it is apparent that the 
classical theory is incapable of accounting for the condi- 
tions associated with the smooth rigid clamping mem- 
ber. It is therefore both desirable and necessary to 
utilize the E. Reissner plate theory? which accounts for 
the effect of transverse shear deformation{ and which 
stipulates the satisfaction of a system of natural edge 
boundary with the physical 
boundary conditions at the edge of the plate. In this 
connection it is pertinent to recall the applications of 
the Reissner theory to plates on an elastic foundation 
carried out by Naghdi and Rowley* and to a class of 
nonlinear axisymmetric plate problems carried out by 


conditions consistent 


one of the present authors.‘ 

In the present paper, with reference to axisymmetric 
problems, by including the effect of the transverse shear 
deformation, expressions are deduced for the stresses 
and displacement (radial rotation) which are valid for 
deflection —-i.e., the region 

At the boundary of the 


the region of constant 

within the smooth clamp. 
clamped region, continuity requirements for the radial 
bending couple and the displacements are satisfied and 
the discontinuity in the shear stress resultant is justi- 
fied by assuming that one of the smooth clamping sur- 
faces exerts a line load on the plate at the edge of the 


clamped region. In addition, an expression for the 


t The Reissner theory accounts for the effect of the transverse 
normal stress, as well as the effect of transverse shear deformation 
In the examples of references 3 and 4, as well as in many other 
examples, the former effect has been demonstrated to be suf 
ficiently small that it may be neglected, at least from a practical 
The same result may be obtained for the example 
of the present paper. Therefore, the effect of the transverse 
normal stress will not be included in the present treatment 


point of view 
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Fic. 1. Loaded plate. 


normal tractions on the plate surface transmitted by 
the surface of the clamping member is obtained. Asa 
specific example, an annular plate with a center annular 
region constrained by smooth clamping surfaces and 
loaded at the outer periphery by a uniform edge load is 
considered. Numerical results afford a means of com- 
parison of the classical theory prediction with the pres- 
ent treatment and indicate the marked significance of 
the present treatment in the case where the radial di- 
mension of the annular portion constrained within the 
clamp is small in comparison with the plate diameter. 


(2) The Basic Equations 


In terms of the dimensionless radius p = r/b, where 
r is the radial coordinate, and 6 is the radius of the 
circular boundary of the clamped region, the basic 
equations of the Reissner plate theory (in polar co- 
ordinates and under the assumptions of axial symmetry) 
aret 


M,' + (M, — M,)/p = bV (pV)' = bpp (1) 
M, = (D/b)(B’ + v(B/p)| My = (D/b)((B/p) + vB’| 
B = —(w’/b) + [12(1 + »)/5Eh|V (2) 


where D = Eh*/12(1 — v?) and prime denotes dif- 
ferentiation with respect to p. In Eqs. (1) and (2), 
M, and M, are, respectively, the radial and circum- 
ferential bending moments; V is the radial shear-stress 
resultant; p is the difference between the normal sur- 
face tractions at the bottom and top surfaces; 6 and w 
are the radial rotation and deflection of the middle 
plane, respectively; / is the thickness of the plate; F is 
Young’s modulus; and vy is Poisson’s ratio. The co- 
ordinate system is chosen so that the deflection is posi- 
tive when measured downward. It also should be re- 
called that the edge boundary conditions consistent 
with the foregoing set of equations are (1) either J/, 
specified or 6 specified, and (2) either V specified or w 
specified. 

Tt In these equations the effect of the transverse normal stress 


has not been included. 
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(3) The Region Within the Clamp 


In the region of the plate constrained within the 
smooth clamp, the deflection, w, is constant. It then 
follows, through appropriate substitution into the first 


of Eqs. (1), that the shear stress resultant, V’, is gov- 
erned by the modified Bessel equation 
V" + (1/p)V’ — [k2 + (1/p?2) ]V = 0 (3) 





where k? = 5(1 — v)(b/h)? 
It then readily follows that 
Ve = Al(kp) + BRyi(kp) 
B = [12(1+ )/S5EA)|[AT (kp) + BK, (Rp) | 


h, k | 
M, waar, A 1)(Rp) — I(Rp) 
5b \ l—p p 


A 


p 
h* | kv l 
i= — fA I\(kp) + - li(Rp) 


5b Li —» 
| kv " . | 
—B Ko(kp) — — Kyi(Rp) | ¢ 
l_—vyp p f 
where J,(kp) and K,(kp) are modified Bessel functions 
of order n. The arbitrary constants of integration, 4 
and 8B, are determined by the specification of either 
M, or B at the circular edges of the region under con- 
sideration, while the second of the edge boundary con- 
ditions is satisfied since the deflection is prescribed. 
The second of the equilibrium conditions, Eqs. (1), 
gives an expression for the difference between the nor- 
mal tractions at the top and bottom surfaces 


k l 
— al Ko(kp) + K(k) | (4) 
— 


— 


p = k/b[Al(kp) — BKo(kp) | (5) 


Let us denote the normal tractions at the top and bot- 
tom surfaces by p; and pf», respectively. Then, in view 
of the assumption that the transverse normal compo- 
nent of the displacement, w, is independent of the thick- 
ness coordinate and the assumed perfect rigidity of the 
smooth clamping surfaces, it may be concluded that 
either p; = 0 po = p, fir = —p po = V, or pi = fe = V 
according to whether p > 0, p< 0, orp =O. In addi- 
tion, it will be noted from the first of Eqs. (4) that in 
general the shear stress resultant, |’, does not vanish in 
the region under consideration. Thus the theory pre- 
dicts that both at a plate edge within the clamp and at 
the boundary of the clamped portion one of the clamp- 
ing surfaces will transmit a line load to the plate. It 
should be pointed out that such line loads account for 
the discontinuity in the shear stress resultant and are 
consistent with the satisfaction of equilibrium. 


(4) An Example 


As an example, let us consider an annular plate of 
inner radius d and outer radius a, concentrically con- 
strained between two smooth rigid circular disks of ra- 
dius 6 (and situated a distance / apart), and loaded by a 
uniform line load of intensity P per unit length at its 
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outer periphery, as shown in Fig. 1. It is convenient a/b. For the regionn < p < 1 the solution is given by 
to introduce the dimensionless radii n = d/b and m = Eqs. (4) and (5). For the region | < p < m the solu- 
tion to Eqs. (1) and (2) may be written 


8 = b’mP/2D[p log p + Cip + (Ci/p)| 
— ee es ae 7 Gi) _ ( 4 c:) - | ‘ ah?P nai 

2D 4 2 2 i (1 — »y)D ; 
YY =maP/o 7 


M, bmP/2((1 + v)(log p + Ci) + 1 — (1 — v)Co/p?)] 

M, = bmP/2[(1 + v)(log p + CG) +r» + (1 — v)(Ci/p?) | 
where C, and C, are arbitrary constants. The two additional constants of integration, arising in the solution of Eqs. 
(1) and (2) have been evaluated from the boundary condition |’,_,, = P and the continuity condition w, -, 0) 
The arbitrary constants, .1, B, C,, and Cy, are determined from the continuity requirements on .\/, and 6 at p x 


together with the boundary conditions 


It is convenient to introduce the functions 


fO = [0/1 -— v0 — hO ( 


> . j (S) 
gd = [0/1 — v|JKoQ) + K,()I 
in which case 
B = [f(kn)/g(kn) |A 
5 2 (1 — p) l— yp C2 l f(Rk . 
A= (¢) : | ( + m) _— logm — EG +- = Ki(b) | 
2\h m l+p m* l+p g(kn) 
, l— vp C2 l 
Ci = - — log m — 
1+ vm? l+ yp (9) 
C: j f(kn) ] (1 + »v) f(kn) ) (; + 
= 4 log f(k) — - g(k) | — I\(k K(k) jlog 2 
m* (; + py T log mi g(kn) g rT (l — “| ie) g(kn) NS) 6 7 (\il — p sis “ 
f(kn} Bs f(kn) _. ) ; 
f(k) — g(k) | + (m? — 1)] i(k) + K(k) |e 
g(kn) g(kn) f 


It should be noted that the clamping member exerts two line loads on the plate: one at p = n of intensity AJ,(kn) + 
BK,(kn) and the other at p = 1 of intensity mP — AJ,(k) — BK,(R). 

It should be mentioned that for the present example (where there is no surface load and the shear stress resultant is 
statically determinant) the predictions of the Reissner theory, employing the usual clamped edge conditions —viz., 
w,-1 = B,=, = 0, are identical to those of the classical treatment insofar as the stresses are concerned. The only 
difference in the prediction of the deflection is contained in the last term of the right-hand side of the second of Eqs 
(6). 

Two limiting cases of the parameter 7 are of interest. As m approaches zero, corresponding to the case of a solid 


plate, the expressions given in Eqs. (9) approach the following values: 


B=0 
| G) (1-—v) P (| =p 4 “ l | 
A = m? — log m — 
2\h m N(R) \1l + vy m* "8 1+ yp | 
l1—vG l (10) 
C, = . — log m — 
1+ vm- l+p 


a l+y js | . 
‘ + log m)f(k) — I,\(k)log m + m*) f(k) + (m? — 1)/(k) 
m- 1+ pv ; l— p 1+ p 


As n approaches unity, in Eqs. (9) the forim of C; is 
unaltered and the expression for C; takes the form ment approach those of the Reissner theory solution 
with boundary conditions appropriate to a simply sup- 


C, = —[(1 + v)log m/(1 — v)(m? — 1)] ; 
ported inner edge. 
Substitution of these values for C; and C, into Eqs. (6) For the case of vy = 1/5 the ratios of w to w, evaluated 
gives a result which is identical to the case where the at p = mand J, to M,, evaluated at p = 1 (w, and 
inner edge of the plate is simply supported. Thus as x M,, being the deflection and radial bending moment, 
approaches unity, the predictions of the present treat- respectively, predicted by the classical theory) as func- 
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Fic. 2. Comparison of deflection with classical-theory 


prediction. 


2 and h/a = 1/5, 1/10, 1/15 are 


The convergence to the case of 


tions of m for m 
shown in Figs. 2 and 3. 
the simply supported inner edge as ” approaches unity 
will be noted from these figures. It will also be ob- 
served that while, regarding the stresses, the improve- 
ment of the present treatment is appreciable only for 
large values of n, regarding the deflection, the improve- 
ment is substantial over the entire range of m (and is 
very substantial as ” approaches unity) and that the 
latter is true even in the case of a very thin plate. 


(5) Conclusions 


It is apparent that in the analysis of plates, portions 
of which are constrained between two smooth, rigid, 
plane surfaces, the utilization of the E. Reissner plate 
theory is essential. The limitations of the classical 
treatment requirements, that both the stresses and dis- 
placements vanish in the region constrained within the 
clamp, are overcome by employing a plate theory which 
includes the effect of 
By accounting for the latter effect, it is possible to de- 


transverse shear deformation. 
duce expressions for the stresses and displacements 
(rotations) valid for the region contained within the 
clamp and to apply appropriate continuity conditions 
at the edge of that region. (It may be observed that, 
at least in the case of thin plates, the use of the classical 
theory for the region not contained within the clamp 
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Satellite Motions 


(Continued from page 802) 
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Fic. 2. 


may give a result which is adequate from a practical 
point of view.) 

On the basis of the numerical results for the axisym- 
ietric example considered, it is clear that, regarding 
the deflection, the improvement of the present treat- 
ment over the classical theory prediction is substantial 
even in the case of an extremely thin plate. It may also 
be concluded that the improvement of the present 
treatment is markedly influenced by the ratio of the 
radial dimension of the clamped region to the radius of 
the plate and that, if this ratio is small (corresponding to 
large values of 7), the improvement is very substantial, 
as regards both the stresses and the deflection. 
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A Further Note on the Laminar Jet Mixing 
of an Electrically Conducting Fluid 
in the Presence of a Magnetic Field* 


K. Toba 

Department of Aeronautical Engineering 
Rensselaer Polytechnic Institute, Troy, N.Y 
November 18, 1960 


Fy rsecges of the laminar jet mixing of an electrically con 
ducting fluid in the presence of a magnetic field has been 
presented by the present writer The magnetic Reynolds num 


ber was assumed to be small. In the following a further dis 
cussion is given for the case when the magnetic Reynolds number 
islarge.t Notation is the same as in ref. 1 

Suppose the magnetic Reynolds number is large such that «€ is 
greater than or close to 1. (The flow Reynolds number is as 
sumed sufficiently large. ) 

For this case we do not assume the applied magnetic field a 
shall seek a field which 


conditions 


priori as a function of x and/or y, but 


vields a similar solution under some prescribed 
Since the magnetic field will be frozen into the fluid particles far 
away, a magnetic boundary layer may be assumed in the jet 
wake region. Thus, V24 in the fundamental equations of Ref. 1 
may be replaced by 324 /dy*®. This is equivalent to the asympto- 
tic expansion employed in a problem of the laminar boundary 
layer of an electrically conducting fluid over a simi-infinite flat 
plate.? In ref. 2 this procedure has been carried out for € varying 

* This research was supported by the United States Ajr Force through the 
AFOSR of the ARDC AF49(638)-23 
in whole or in part is permitted for any purpose of the United States Govern 


under Contract No Reproduction 


ment 
including the discussion of the “‘plug 
K. Kusukawa informed the author of 


in which the same problem is treated together with the axisymmetric 


t+ After the present formulation 
ging’ was done, Dr 
Ref. 3 


case 


phenomenon 


The solution is found analytically by matching the outer and inner 


solutions 
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1.0 with the boundary conditions: P 
| ——— Fi 
Ph = f=0, F=G=0) o| 
Gi ys \ \\\———— +4 t— o, f= G’ =0} } i 
| ee — ° 
Poe — .95 : ? , : ? 
8 1G Numerical integration has been carried out by means of an (I 
H . 1 —__|__¢ analog computer for various values of 8 and e land 4.§ The y 
ay Fh ER results are shown in Figs. 1 and 2. F’ and G’ are to be recog “a 
= rs : - —— nized as u/Ugxis and H,/H, axis, respectively. In Ref. 3 numeri 
rome: cal results are presented in curves of not only u and H, but also 
the streamlines, magnetic-field lines, and the current density 
with e l 
Comparison of the present results with those of Ref. 3 shows a 
4 - —— marked difference in the u and H, profiles. According to Ref. 3 P 
16 the curves approach zero more slowly as # increases; whereas 2 
the present results show that there is no such remarkable differ 01 
ence, regardless of the values of 8. This discrepancy warrants 
further investigation. H 
The possible numerical error may be estimated by the relation Ai 
( F)? B(G)? + 1 7) H 
Ja 
which is derived from Eq. (5). The quantities in parentheses 
denote the corresponding asymptotic values. The present com- 
1 j___. 1. E putation gives errors roughly less than 10 percent for all cases 7 
8 16 ee 
plotted except for 8 = 0.95 W 
Lastly, it should be pointed out that according to the machine hes 
computation the solution seems to exist even for B > | Since th 
the error is increasing close to 8 1, it is admittedly difficult to 
~ « . e . - is 
draw any definite conclusion from the numerical computation; 
it still appears that the proot for the existence of the ‘“‘plugging” 
e ° o ae 
phenomenon needs more perfection, except fore = © re 
f 
REFERENCES ‘i 
= ——e ———as 16 gE 1 Toba, K., Laminar Jet Mixing of an Electrically Conducting Fluid in the 
- . , - ee wh 
Fig. 1. Numerical results fore = 1. F’ denotes u/taxis and G 
denotes H,/Haxis § Thanks are due to Mr. W. McClintock, without whose cooperation this iro 
work would have been delayed unduly a; 
an 
tel 
col 
from 0 to ©. However, as seen from Eq. (3) of Ref. 1, for very 
F ‘ . : om : oe ae 
small values of e diffusion may dominate over convection, so that B=O wen 
L nom & { 
we may not be justified in neglecting 0?4/0x? compared to —/f ~ spt 
0°A/dy? as for small magnetic Reynolds number. This is not 4 ~ ing 
e . . . . - > 
so for the present case. Hence, if a similar solution tie 
. EE . is 4 A 
= £1 3f(n)! ‘ it ° (ra 
v hy - 9 = y/Bx?3 (1) it at 
A = xV3g(n) | Bs ea slo 
: ; P , P : , eh of 
exists, it must satisfy the following differential equations: 
wa 
Fig a h ” am f' 2 B( gg” ae g! 2) ] a 
la = fg’ = od 9 it 1 = i = 
De eae: (1/e) g ) (2) 8 16. € req 
with the boundary conditions I 
- : 7 ) wit 
= 0 =f" = g"=0 mes 
7 , J Jj ‘ g » (3) can 
7 ©, 7 =2 = 0 \ : z 
' of 1 
js . 6 .5r 
It has been shown? that for a laminar boundary layer over a the 
semi-infinite plate of an electrically conducting fluid under the e% que 
parallel uniform magnetic field, the flow is ‘‘plugged”’ as B —> 1, te of § 
e - e » ; a = 
or the velocity field and the magnetic field are annulled at this cree 
critical value of 8 and no useful solution can be found for 8 > 1.f O en | ni iner 
For the present problem the “‘plugging’’ seems to occur as B > 1; 1.0 8 16 & criti 
this can be most easily checked for ¢ ~ © by noticing the can- I 
cellation of the inertia terms by the magnetic body-force terms. crea 
Eqs. (1) and (2) can be further simplified by letting quir 
ta . an 
n = t/a, f(n) = 2aF(E), 2(n) = 2aG(é) (4) ‘ ; 
tan 
where a is an arbitrary constant. We then have lr 
F'+ F = pG?+1 . only 
ye - ~ (o an 
F'G — FG' = (1/2€)G’” ) 8 twic 
; J 1 5/2 
<7 
_ 16 & Iner 
t Although not explicitly mentioned,’the same phenomenon may be ae 
observed for a wedge flow.‘ Fig. 2. Numerical results for e } ' 
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Parametric Study of the Influence of 
Propellant Sloshing on the Stability 
of Spacecraft 


Helmut F. Bauer 

Aeroballistics Division, George C. Marshall Space Flight Center, 
Huntsville, Ala. 

January 10, 1961 


a DANGER of unstable modes resulting from propellant 
sloshing for rigid or nonrigid spacecraft is well known. 
While the exact treatment of these phenomena, especially in the 
latter case, leads to rather unwieldy equations which do not lend 
themselves to easy generalization, very useful general results can 
be obtained from a simplified stability investigation. 

The control equation in this study is idealized by neglecting 
derivatives of 8 which produce increasing phase lags with increas- 
ing frequency, but are not important for the basic argument. It 


1s 
B = angi + agi + {boa; or/and goA;} 


where B is the gimbal angle, ¢; is the indicated attitude deviation 
from a space-fixed reference (as measured from a position gyro), 
a; is the indicated angle of attack, A; is the indicated acceleration, 
and do, a1, bo, and ge represent the gain values of the control sys- 
tem. (Note: half of the thrust was considered available for 
control purpose. ) 

The forces and moments due to propellant sloshing were 
described by an equivalent mechanical model, which essentially 
consists of a nonsloshing mass with a moment of inertia and 
spring-dashpot-mass systems, each of which represents one slosh- 
ing mode. For simplification only one tank and one sloshing 
mode were considered and the influence of the various parameters 
(ratio of sloshing propellant mass to total vehicle mass, natural 
slosh frequency, control frequency, control damping, gain factors 
of the attitude-control system, accelerometer location and gain) 
was investigated. The stability boundary was determined in 
terms of the amount of damping of the propellant in the tank 
required for various tank locations (slosh-mass location). 

For a rigid spacecraft with attitude control only (position gyro 
with differentiating network), the danger zone where instability 
can occur if not properly baffled is between the center of gravity 
of the spacecraft and the center of instantaneous rotation. If 
the slosh frequency is below the control frequency (natural fre- 
quency of pitch mode) the danger zone increases aft of the center 
of gravity and stronger baffling is required for stability. In- 
creasing the control damping in the subcritical region leads to an 
increase of required baffles and to a decrease of baffles in the super- 
critical region. 

Decreasing the gain value da) at constant control frequency in- 
creases the danger zone slightly to the aft of the vehicle and re- 
quires stronger baffling. Furthermore a decrease in baffling 
can be obtained by increasing the slosh frequency (change of 
tank form). 

In the treatment of an elastic spacecraft with attitude control 
only, the results indicate an enlarged danger zone and about 
twice as much required damping for stability (slosh frequency 
5/2 times and bending frequency 5 times the control frequency). 
Increasing the control frequency enlarges the danger zone 
toward the tail of the vehicle and requires considerably more 
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Fic. 1. Stability boundary 


damping in the tanks. Increase of the control damping has the 
same effect as mentioned above, except that about two to three 
times more baffling is needed due to bending. Increasing 
bending frequency results in a decrease of the danger zone and 
amount of damping and approaches the values of a rigid space 
vehicle. The change of parameters previously mentioned 
shows the same trend as in the rigid space vehicle case with the 
exception of increased baffling. 

Additional artificial stabilization through the control system, 
such as angle-of-attack and/or accelerometer control, can en 
hance the stability behavior if the gain values and sensor char 
acteristics are properly chosen. Additional angle-of-attack 
meter control in a rigid space vehicle indicates the same danger 
zone as in the rigid space vehicle (with gyro and differentiating 
network only) but reduces the required baffling for stability by 
The trend of the other parameters is 
Additional 
accelerometer control (perpendicular to the axis of the vehicle 
alleviates the influence of propellant sloshing considerably, if the 
location of the accelerometer, its characteristics, and gain value 
are properly chosen. For a gain value ge = 1/g (g being the 
longitudinal acceleration of the spacecraft) the danger zone is 
reduced to a very small region about the center of instantaneous 
rotation, where very little damping is needed for stability. The 
increase of the control frequency results in an increase of the 
danger zone toward the aft of the vehicle and an increase of re- 
quired damping. The required baffling, however, is far below 
Control damping 


a factor of two to three 
very similar to the previously mentioned cases 


those values without accelerometer control 
increase again requires an increase in baffling in the subcritical 
region, while it requires a decrease in baffling in the supercritical 
region. For slosh frequencies below the control frequency the 
danger zone increases aft of the center of instantaneous rotation 
The location of the accelerometer for a proper gain value g. = 
1/g is of minor importance, but has definitely a strong influence 
for different gain values. If ge = 2/g, the location of the ac- 
celerometer aft of the center of gravity would make the vehicle 
unstable for all tanks behind the center of instantaneous rotation; 
while for an accelerometer location in front of the center of 
gravity, all tanks in front of the center of instantaneous rotation 
would make the vehicle unstable, or would require very strong 
baffling for stability. For values go < 1/g the danger zone lies 
between the center of gravity and the center of instantaneous 
rotation, while values g2 > 1/g shift the danger zone forward of 
the center of instantaneous rotation (for accelerometer location 
slightly forward of the center of gravity of the vehicle). These 
results are valid for an ideal accelerometer or real accelerometer 
with large natural frequency. For decreasing natural frequency 
of the accelerometer the stability boundary shifts above the x,- 
axis again, thus demanding baffling in the tanks. 

In Fig. 1 the stability boundary is represented for the four 
cases mentioned, the abscissa being the location of the sloshing 
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mass and the ordinate being the required damping in the tank. 
It is assumed that 16 percent of the spacecraft’s mass is sloshing 
at the location of the slosh mass. Values y, above the curve 
indicate the stable region. It can be seen that bending increases 
the danger zone and requires larger damping of the propellant 
Furthermore, the effect of angle-of-attack meter and accelerom- 
eter control with go = 1/g is clearly expressed in the respective 
curves. 

The potential hazard in control due to propellant sloshing in 
space vehicles can be eliminated by proper choice of tank form 
(slosh mass ratio decrease and slosh frequency increase), proper 
selection of type, location, and gain values of control sensors, 


and—as a last resort—by baffles. 


Critical Frequencies in the Stagnation Region 
of a Shock Layer; 


R. J. Leite and N. £. Hawk* 

Radiation Laboratory, Department of Electrical Engineering, 
University of Michigan, Ann Arbor, Mich. 

January 9, 1961 


7 FACILITATE the estimating of stagnation-region plasma and 
electron, neutral-particle collision frequencies, a plot of alti- 
tude versus velocity with equilibrium plasma and collision fre 
quencies indicated parametrically was developed. Since the gas 
in the stagnation region is relatively dense and only weakly 
ionized, a so-called critical frequency is introduced to account for 
the influence of these elastic collisions. Alteration of the plasma 
frequency, when collisions are considered, is indicated on the 
above-mentioned plot. 

The term ‘“‘plasma,’’ as used here, refers to an electrically 
neutral, ionized gas consisting of equal populations of electrons 
and ions, and neutral particles also. The assumption of elec- 
trical neutrality implies that the effects of space-charge forces, 
namely internally generated electric fields, are not important at 
distances comparable to the shortest length of interest. In the 
case of electromagnetic radiation, this length is the wavelength 
of the radiation considered. It is well known that a character- 
istic frequency, called the plasma frequency, w,, is associated 
with each plasma, namely, 


» 


wp? = ne*/me (1) 


where » is the electron number density, ¢ is the charge of an elec- 
tron, m is its mass and e is the dielectric constant of free space 
Physically, wp can be thought of as the frequency at which elec- 
trons will oscillate about their equilibrium positions after being 
displaced by an external force. Upon being displaced, the elec- 
trons are subjected to a restoring force due to the electric field 
between them and the ions, which move very little because of their 
greater mass. This restoring force will tend to return the elec- 
trons to their equilibrium positions and will cause them to oscil- 
late, in an undamped manner, about these positions, in the ab- 
sence of a magnetic field. It can be shown that the influence of 
free electrons is always to reduce the dielectric constant below 
the free-space value, namely, 

e’ = e[1 — (w)/w)?] (2) 


where e’ is the effective dielectric constant of the plasma and w 
is the frequency of an incident electromagnetic wave. 
Determination of electron density in the stagnation region of 


+ This research was supported by the Advanced Research Projects Agency 
and the U.S. Army Signal Research and Development Agency under Con 
tract No. DA 36-039 SC-75041, ARPA Order No. 120-60 Project Code No 
7700 

* Research Engineer and Assistant in Research, respectively 
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the shock layer during re-entry requires reference to some model 
atmosphere.! In addition, real-gas effects must be considered 
in the evaluation of properties? and electronic states® of the air 
in that region. In Ref. 2 the Rankine-Hugoniot shock-wave 
equations are solved with dissociation and ionization included 
Using these shock-wave results and the data in Ref. 3, the equi- 
librium electron density as a function of velocity and altitude cay 
be computed easily. Using the relationship 


fo = 8.97 Vn 3 


where fp = w,/2z is the plasma frequency in cps and x is the 
number of electrons per cubic meter, loci of constant plasma fre 
quency are drawn on the altitude-velocity plot shown in Fig. | 
Superposition of any altitude-velocity trajectory on this plot will 
indicate the plasma frequency in the stagnation region at points 
along the path of the re-entry body. 

In general, a plasma will reflect electromagnetic energy which 
is radiated at frequencies less than f, and it will absorb this energy 
when the radiating frequency is equal to or greater than f,. At 
frequencies only slightly greater than f,, attenuation of the inci- 
dent electromagnetic wave is quite large; however, as the propa- 
gation frequency becomes increasingly greater than f, the effect 
of the free electrons becomes less pronounced and conditions 
approach those of free space. This can be readily seen by re- 
ferring to Eq. (2), which shows that ¢€ is approached as propaga 
tion frequencies become greate> than f, 

The above would be true in a collision-free gas 
layer however, the gas density is of such magnitude that free- 


In the shock 


electron, neutral-particle collisions must be considered during at 
least part of the re-entry flight. Collisions between free elec- 
trons and positive ions can be disregarded, however, because 
of the degree of ionization is usually small—of the order of several 
percent. Hence, the above discussion of the propagation of elec- 
tromagnetic waves in the stagnation region must be modified 
when the electron, neutral-particle collision frequency, v, begins 
to approach f,. When this happens, the classical definition of 
plasma frequency [Eq. (1)] will remain unchanged, although 
its physical significance becomes doubtful. In other words the 
ordered motion of the electrons tends to be made random by the 
clectron, neutral-particle collision processes and this, in turn, in- 
fluences the electromagnetic properties of the plasma 

Referring to Whitmer,’ we see that ¢’ of a plasma in which 
collisions are considered is 

Wp" _ Wp*(v/w) 


c = 4 1 — " i? s (4) 


9 
le lh ad rr er 


where y is collisions per second. The inclusion of v has resulted 
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(2 x fp)? = (ne?/me); 


vy = 3.6 X 10" po/po (72/1,000)"2, f.2 = fp? [1 — (v/wp)? 


Typical trajectories® of several types of missiles are also shown 
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jn a complex value of e’, thereby indicating that the plasma now 
has some additional loss properties which further attenuate any 
propagating electromagnetic wave. In addition, it can be shown 
that the electron velocity is no longer 90° out of phase with the 
electric field of the incident electromagnetic wave as it was in the 
collision-free case. This means that the free electrons can now 
absorb energy from the incident electric field, thereby increasing 
the degree of ionization of the gas because of the occurrence of 
higher-energy interactions during electron, neutral-particle 
collisions. 

Eq. (2) shows that in the collision-free case the critical fre- 
quency of the gas (that frequency at which absorption is maxi- 
mum) is equal to w,. When collisions are considered, Eq. (4) 
indicates that the critical frequency of the gas is dependent upon 
both w, and »v. This is seen readily by setting the real part of 
Eq. (4) equal to zero, namely, 


we? = w,*[1 — (v/w,)?] (5) 


where w, is the critical frequency of the gas with collisions. We 
see that when y is small compared with w,, Eq. (5) vields results 
for the collision-free case 

Eq. (5) has been evaluated in the shock layer and the results 
are shown in Fig. 1, where f. = w-/27. Collision frequency was 
determined from the relationship 


vy = 3.6 X 10!" po/po( T2/1,000)! 2 (6) 


where p2 and 7» are stagnation-region values of density and tem- 
perature, respectively, and py» is standard sea-level density. It 
should be noted that when + is less than 10 percent of f,, the propa- 
gation conditions of a collision-free gas may be assumed for all 
practical purposes. On the other hand for greater values of 
v, the f. of the gas decreases with increasing v, thereby permitting 
additional frequencies, namely wp > w > w:, to propagate through 
the gas. In the limit, when v = wp», at least in principle, all fre- 
quencies propagate through the gas with some degree of attenua- 


tion 


REFERENCES 
Minzner, R. A., Champion, K. S. W., and Pond, H. L., The ARD¢ 
Vodel Atmos phere, 1959, Geophysics Research Directorate, AFCRC, ARDC, 
Air Force Surveys in Geophysics, No. 115, AFCRC-TC-TC-59-267, August 
1959 
2 Hochstim, A. R., Gas Properties Behind Shocks at Hypersonic Velocities, 
I. Normal Shocks in Air, Convair, San Diego, Physics Group, Rep. No 
ZPh(GP)-002, Jan. 30, 1957 
3 Gilmore, F. R., Equilibrium Composition and Thermodynamic Properties 
of Air to 24,000°K, The RAND Corporation, RM-1543, Aug. 24, 1955 
Whitmer, R. F., Principles of Microwave Interactions with lonized Media 
The Microwave Journal, Vol. 2, No. 2, pp. 17-19, Feb. 1959 
®’ Rose, P. H., and Stark, W. I., Stagnation Point Heai-Transfer Measure 
» 


ments in Dissociated Air, Journal of the Aeronautical Sciences, Vol. 25, No. 2, 


pp. 86-97, February 1958 


Analytical Expression for the Roll Rate of 
Aileron-Equipped Ballistic Re-Entry Bodies 


Martin R. Fink 

Supervisor, Missile Aerodynamics Group, 

United Aircraft Corporation Research Laboratories 
East Hartford, Conn. 

March 3, 1961 


List OF SYMBOLS 


A drag reference area, ft? 

b reference length, ft 

Cp = drag coefficient 

C.(6) = rolling-moment coefficient caused by ailerons 
Clip = damping-in-roll coefficient, OC) 0(wb/2V) 

Tez = moment of inertia in roll, slug-ft? 

ko = drag parameter defined by Eq. (5) 

K = damping parameter defined by Eq. (7) 

m mass of re-entry vehicle, slugs 

Ss = rolling-moment reference area, ft? 


t = time, sec 
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V = velocity, ft. sec 

VE = velocity at start of re-entry, ft, sec 

¥ = altitude above surface of the earth, ft 

8 = exponential altitude parameter defined by Eq. (3 22,000 ft -! 

OE = angle between local horizont®l and re-entry velocity, positive 
downward 

p = atmospheric density, slug-ft* 

po nominal low-altitude density defined by Eq 0 0034 slug 
ft# 

w roll rate, rad ‘sec 

We equilibrium roll rate, rad /sec 


— re-entry bodies may require an angular velocity 
about the longitudinal axis to minimize dispersion caused by 
inadequate static stability and nonzero trimmed lift coefficient 
A rolling velocity can be imparted either by aileron surfaces 
affixed to the body, which produce a roll rate which increases 
during re-entry, or by reaction jets, which may be heavier and 
less reliable than fixed aileron surfaces but which can produce 
the desired roll rate prior to re-entry. In order to evaluate these 
two methods, an analytical expression is required for the variation 
of roll rate with altitude for ballistic re-entry bodies which utilize 
aileron-induced roll. Such a relation is derived herein from 
analytical studies of ballistic re-entry trajectories 

The angular acceleration in roll, dw/dt, is determined by the 
equation 


I,1(dw dt) = C(6)(p/2) V2Sb 4 C) (p/2 V2Sb(wb/2V 1) 


where C;(6) is the rolling-moment coefficient due to control-panel 
deflection and C;, is the damping-in-roll coefficient. The tra 
jectory of a typical nonlifting ballistic re-entry body is essen 
tially linear at the altitudes at which roll is developed, and the 
re-entry trajectory analysis given in Ref. 1 can be utilized 
Thus, the rate of change of altitude with time can be approxi 
mated by 

— V sin 6; 2) 


dy/dt 
and the variation of atmospheric density with altitude can be 


approximated by 


p put ] 
From Eq. (13) of Ref. 1 
! By 
V = Vee ? ' 
where k, = (po/B sin 6¢)(CpA/m 5 
Eq. (1) can then be rearranged as 
ds : Sb koVe m ke — By 
+ Ke~?¥ wo = Clb —e~h¥- 2 6 
d(—By ~~ 2 Cea 
where 
K= - Ci, (.Sb?/1,2)(po/4B8 sin 6, (7 


The aerodynamic coefficients are assumed constant and the roll 
rate is assumed equal to zero prior to re-entry Eq. (6) can then 


be integrated to give 


“(8 Sb y m ky ( — oye? ") 
» = Ci - Sa er am € 2 8) 
nie al ae ae 


in rad/sec, which reduces to 


w = C,(6)(Sb/I,2) Ve(m/CpA (, coe ) 


if the damping-in-roll coefficient, Ci, is assumed equal to zero 
The calculated variation of roll rate with altitude is tabulated 
below for a representative ICBM re-entry body which has an 
average roll control moment per unit dynamic pressure, C;(6).Sd, 
of 1 in.*. This value of roll control moment includes a large 
loss of control effectiveness caused by immersion of the panels 
within the boundary layer. The tabulated values of roll are 
based on a moment of inertia, J,,, of 15 slug-ft®, a ballistic coeffi- 
cient, mg/CpA, of 1,000 psf, a minimum-energy 5,500 nm tra- 
jectory, and damping in roll such that A/& is equal to zero (no 


damping in roll) and 0.1. For this case, the roll rate increases 
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TABLE 1 
Altitude, Roll Rate, o, rpm 
y, ft K/ko = 0 K/ko = 0.1 
200,000 0.1 0.08 
150,000 1.0 0.7 
100,000 9 6 
50,000 79 51 
0 270 114 


by an order of magnitude with each 50,000 ft decrease of altitude 


except at the lowest altitudes. 

The equilibrium roll rate, at which the rate of change of roll 
rate with time would be equal to zero, can be obtained from 
Eq. (1) 

Ws = —[C1(6)/Ci,] (2 V/b) (10) 


By substitution of Eqs. (4), (5), and (7) into Eq. (10), the equili- 
brium roll rate is given by 
ko —By 


we = Cy(6)(Sb/Tzr)(m CpA)(ko/2K)Vre_ i (11) 


and the ratio of instantaneous roll rate to equilibrium roll rate is 


ko -\ —8y 
- K )e , ‘ 
W/We = [2K /(Ro + 2K)] [ .G ) alia J (12) 


Roll rates calculated by use of Eq. (8) can be utilized for deter- 
mination of aileron rolling-moment coefficients required to pro- 
duce a specified roll rate at a specified altitude, for comparison 
with other methods of producing rolling velocities on ballistic 


re-entry bodies. 
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Derivation of Angular-Velocity Ratio 
for a Particular Despin System 
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SYMBOLS 


IJ = moment of inertia of payload or main body about spin axis 

I~ = combined moment of inertia of ali cables in radial position about 
spin axis 

L = length of individual cable 

m combined mass of all small masses 

M = combined mass of cables and small masses 

R = radius about which mass-cable units are wound; also, radius to the 


cable-release point 
= final angular velocity of payload when mass-cable units are in the 


o = 
radial position 

w) = initial angular velocity of payload and mass-cable units prior to 
despin 

Q angular velocity of mass-cable units in radial position about spin 


axis 


INTRODUCTION 

ly RECENT YEARS the use of spin stabilization has been an ade- 

quate solution for many of the stability problems of rocket- 
boosted payloads. However, in some cases the high spin rate 
necessary for stability during the vehicle’s ascent through the 
atmosphere becomes an undesirable parameter once the vehicle 
is essentially out of the atmosphere. In such cases another 
technique short of complete attitude stabilization has been to 
despin the payload to either some lower desired spin rate or 


zero rate. 
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One of the means successfully used to despin a payload has 
been to release small masses attached to cables wound around the 
payload. As the cables unwind and the masses move away from 
the payload, the tension in the cables applies a torque to reduce 
the payload spin rate. Finally, when the mass-cable units are 
released from the payload, there has been an angular-momentum 
transfer to these units resulting in decreased angular velocity of 
the basic payload 

It is the purpose of this note to present a derivation and brief 
discussion of equations relating to the despin results for a system 
such as that described above. Only the relationship between 
initial and final spin rates will be considered; no equations of 
motion or dynamics of the process will be presented. Ref. 1 
concerns the despin problem and other mechanical means for its 


accomplishment. 


DERIVATION 

The following assumptions are made in the derivation of the 
ratio of the final spin rate to the initial spin rate. 

(a) The cables are straight in the radial position of release, the 
end of the despin process. 

(b) The cable lengths are equal and constant. 

(c) The moments of inertia of the small masses about their 
centers of gravity are negligible. 

(d) There are no external forces—e.g., aerodynamic—acting 
during the despin process. 

(e) Angular momentum and energy are conserved. 

The conservation of angular momentum requires that the 
initial angular momentum of the payload and wound mass-cable 
units be equal to the final angular momentum when the mass- 
cable units are in the radial release position: 

(MR? + I)wo = [m(L + R)? + [JQ + Iw (1 

Similarly, the conservation of energy gives 

9 


(MR? + I) (wo?/2) = [m(L + R)? + I,|(Q2/2) + I(w?/2) (2 


Eliminating 2 from Eqs. (1) and (2) and solving for (w/w ) re- 


sults in 
A(w/wo)? + Blw/wo) + C = 0 (3) 
where A = J? + I[m(L + R)? + J,] 
B = —2] (MR? + I) 


C = (MR? + I)[((MR?2 + I) — m(L + R)? — I] 


In the quadratic solution of Eq. (3) for (w/wo) the negative sign 
of the radical is required. 

It is evident that Eq. (3) can be rewritten to solve for the mass 
and/or cable length required to give a specified despin ratio, 
(w/wo ) 

DISCUSSION 

Several other related points of interest are observable from the 
preceding discussion. It is evident that the above equations 
apply for any number of mass-cable units provided they are 
identical. For applications where the mass of the cable is small 
compared to the despin masses, or where less accuracy in the 
final spin rate might be satisfactory, Eq. (3) can be used with 
simpler coefficients obtained by setting 7, = Oand M = m. For 
the special case of despinning to zero final rate, the relationship 
required between parameters is given by C = O from Eq. (3), 
whether the cables are included or not. 

Another result of interest and of some usefulness can be de- 
rived by considering the case where the mass of the cables is 
neglected and the release is made when the cable is tangential 
to the body, rather than radial. If D is the cable length when 
the release is tangential, the rather simple expression for despin 
ratio is 

(w/w) = (I + mR? — mD?)/( + mR? + mD?) (4 


If a comparison is made between the cable lengths required to 
despin to w = 0 from Eqs. (4) and (3) (imodified to omit effects of 
cable mass), it is noted that LZ = D — R. A further approxima- 
tion which can be shown to be quite good here is to use Eq. (4 
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to determine the effective cable length D, then make the actual 
length L = D — Rand release at the radial position. It is em- 
phasized again, however, that these simplifications do neglect 
the mass of the cables, and for many applications this is not 


justifiable 


ILLUSTRATIVE EXAMPLE 

As an example illustrating the effect of cable mass consider a 
payload having a moment of inertia about its spin axis of 2 slug- 
ft?. The despin arrangement is assumed designed such that 
R 1 ftand L = 10 ft, but it is required to determine the mass 
necessary to obtain a despin ratio of 0.1. Two cables are as- 
sumed. If initially it is assumed that the required ratio is a 
maximum desired value, then Eq. (4) might suffice. The mass 
0.0136 slug. If it is further assumed 


required from Eq. (4) is m 
‘ slug/ft, 


that the mass per unit length of the cable is 4.75 (10) 
then Eq. @8) would show the actual despin ratio to be —0.0115 
i.e., the actual final spin rate is much less than the 0.1 limit and 
the direction of the residual spin is opposite to that of the initial 
spin. The large difference noted here is a direct result of the mass 
of the cables, which for this case was about 70 percent of that of 
the despin weiglits. Had it been more important to get exactly 
the despin ratio of 0.1, Eq. (3) should have been solved for m with 
the result m = 0.010 slug 
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Wind-Tunnel Interference for Wing-Body 
Combination 
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Mathematics Dept., Faculty of Science, Alexandria University, 
Alexandria, Egypt 

January 6, 1961 


¢ bw FLOW about a finite wing in free air is largely determined 
by a horseshoe vortex system consisting of a bound vortex 
lying along the quarter-chord span and a pair of vortices trailing 
from points near the tips. When the wing is placed within a 
wind tunnel, an interference is produced at the wing due to the 
limited extent of the flow by the tunnel walls. 

The determination of this interference for a wing alone is now 
well founded! and there is no question as to the validity of the 
obtained results. For many years, however, the same inter- 
ference has been applied to wing-body data, even though this 
practice has never been justified. It is meant here to discuss 
such a practice. For this purpose a wing attached to a circular 
body is placed first inside a circular tunnel and then a rectangular 
one. It will be assumed that the wind stream in the tunnel be- 
fore inserting the model is uniform and the trailing vortices are 


parallel to the tunnel axis 
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We consider first the case when the previously mentioned wing 
body combination is placed inside a circular tunnel of radius R 
Let the strength of the equivalent horseshoe vortex be I’, and its 
span be 2s. Hence, the cross section of the tunnel through the 
wing (Fig. 1) will be at right angles to two semi-infinite vortices 
of strength T, —I° equidistant from the body center and at a 
distance s from it. If the x-axis is taken along the wing with the 
origin at the body center, the flow will be symmetric about the 
y-axis, and the flow on one of its sides may only be considered 

The transformation 


¢ log z/r(—mw < are z < gz, § & + ty l 
where ¢ is the radius of the body, will map the region bounded by 
Zz R, |e r, and the y-axis onto the region bounded by the 
rectangle whose corners are ¢ +in/2, log R/r + iw/2 (Fig. 2 
The inside of the rectangle can then be mapped on the right half 
of a Z-plane (Z = X +7Y) defined by 
P 4 sec (At, R (2) 
where \ and & are such that the corresponding complete elliptic 
integrals A, K’ are given by 
K/2 log R/r K'/x r/2 (3) 
The vortex at s = s will be transformed to a vortex of the same 
strength at Z = sec A log s/r 
The complex potential in the Z-plane is therefore 
ly a s iT ; s 
w= log | Z — sec A log - log | Z + sec A log (4) 
4 r 4a r 


Substituting for Z in terms of z, the induced velocity in the z-plane 


can be obtained ; on the wing this gives 


< 
~~ 
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r dr p 
uu 0 i cos \ log p — cos dA log (5 
trs p pi- 
where p v/S, pi r/s. 


It is evident that when the wing-body combination is placed 
in free air, the velocity induced at the wing is 


r 1 pi" 
ul 0), Y= — 
2s \p? — 1 p> — pi! 


net velocity induced at the wing by the tunnel wall is 


(6) 


Thus the 
v’=v—y (7) 
If lis the flow velocity in the tunnel before placing the model, 
S the wing area, and C,, the wing lift coefficient, then the strength 
of the horseshoe vortex [ = C,(SU/4s) and the interference 
can be put in the form (SC;,/C)é where C 


A 
cos X log = 
p 
p 4 
cos A log a (S 
pi~ 


The values of 6 as represented in Eq. (8) have been plotted in 
0, 0.1, 0.2, 


angle Aa(=v', 
rR? and 


Fig. 3 for values of s/R = 0.3 and 0.7, and of r/s 
and 0.3. 

In considering the case of a rectangular tunnel, let the above 
wing-body combination be placed symmetrically with respect 
to the tunnel center with its span parallel to its breadth (Fig. 4) 
and assume that the radius of the body ¢ is small compared with 
the tunnel breadth 6 and height A. 

Transforming toa ¢-plane given by 


f=2:- (¢=é&+% (9) 


the body will be mapped on the part of the n-axis where —2r < 


n < 2r and the vortices trailing from s = +s to vortices of the 
same strength at ¢ = +) where () = s — r?/s. Since ¢ is small 


* This correction has been obtained by C. B. Smith in the form of an 


9 


infinite series ? 
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compared with 6 and h, the transformation will leave the tunnel 
It is evident that the flow on one side 
This 
flow is now similar to that in the ¢-plane of the previous section 
Therefore, if \ and the moduli k, k’ are such that the correspond- 


walls unaffected (Fig. 5 
of the y-axis, the right one say, need only be considered. 


ing complete elliptic integrals A, A’ are related by 


K/b K'/h = X/2 (10 


the induced velocity in the ¢-plane can be obtained in the form 


hl 


—u + A[cos A(f — fo) — cos ACE + Fo)] (11 
On the wing in the z-plane it gives 
u = 0, v PA 47} cos AS(p — 1){1 + (p;? p } 


cos As(p + 1)[1 — (p1?/p)] 5 [1 + (p12/pe (12 


The interference angle Aa can be calculated as before; it can 
be put in the form (SC;,/C)6 where C is the tunnel cross-sectional 
area and 


bh 2 2p? 
= - | = = AS = As(1 p (, + 
1 ( \ 


167s? — p* p? — p;' 


=) cos As(1 + p) (: _ yt (: re *) (13 
p p f p* 


When s is small, this expression can be approximated to 


6 = (KK'/12)(2 — k2)(1 — p,?)[1 + (p12/p2)] (14 


The distribution of 6 along the wing span for several typical 
shown in Fig. 6. 

In Eq. (14) it can be proved that KA’(2 — k? 
value whether A/K’ b/h or 2h/b. 
rection [Eq. (14)] is the same for a tunnel of height / and breadth 


has the same 


This means that the cor- 


2 and breadth 1/2h; and it can be 
at b/k = 2. This 
|, a quantity 


height 5/, 
6 has a minimum 
minimum value equals 0.119(1 
whose value at any point of the wing span is slightly less than its 


b as for one of 
deduced that value 
— pi®)[1 + (o1?/p? 
corresponding value were the tunnel circular of the same cross- 
sectional area. It seems, therefore, that for small wing spans, 
the best working tunnel of given cross-sectional area is rectangular 


of breadth 4/2 its height. 
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It is evident from the whole set of plotted graphs that the 
interference for wing-body combinations due to the tunnel walls 
is significantly different from the wing alone and that the differ- 
ence increases with p; In particular, this interference is in all 
cases greatest at the wing root, thus tending to produce a root 
stall rather than a tip stall, contrary to that produced by the wing 
alone 

It is felt that these results are particularly significant insofar 
as they affect the interpretation of stall studies. The results 
in general show that the wing-alone corrections should not be 
applied to wing-body data for any purpose, particularly for large 


values of p; 


REFERENCES 
Glauert, H., The Elements of Aerofoils and Airscrew Theory, Cambridge 
University Press, 1948 
Smith, C. B., Wind-Tunnel Wall Corrections for Wing- Body Combinations, 
Journal of the Aeronautical Sciences, Vol. 16, No. 4, pp. 237-242, April 1949 


A Simple Method for the Analog Computation 
of the Mean-Square Response of Airplanes 
to Atmospheric Turbulence 


B. Etkin 
Institute of Aerophysics, University of Toronto 
February 20, 1961 


SYMBOLS 


= wing mean chord 


} reduced frequency 

I turbulence scale 

u reference flight speed 

Wg (t) downwash component of atmospheric turbulence 
ag (t) [—wegy(t) ]/u 

o = intensity of turbulence (rms velocity) 

w = circular frequency 

Qi wave number 


I MAY BE the case when dealing with the response of an airplane 
to random atmospheric turbulence that the mean-square re 
sponse is the only information desired. If so, a method of analy- 
sis which requires first the calculation of the frequency-response 
functions and the output power spectra! is wasteful of effort 
Analog computation offers a means of bypassing these steps, 
and proceeding directly to the desired mean-square value. One 
method of doing this, which utilizes random signals for the input, 
has been described by Mazelsky and Amey.? It is applicable to 
nonlinear as well as to linear systems, and can employ non- 
Gaussian inputs 

For systems which are linear, the methods of analog computa- 
tion given by Laning and Battin*® can be adapted to the gust 
response problem. These do not require the generation and 
filtering of random signals, as does Ref. 2, but still involve some 
complications and practical difficulties 

The method developed herein, applicable to linear systems, is 
considered to be simpler to apply than any of the above. It 
hinges on the fact that the amplitude spectrumf of the simple 
exponential function Ae 7 gives a good representation of the 
power spectral density of atmospheric turbulence. Hence the 
response of a linear system to this easily-generated exponential 
input can be used to calculate its mean-square response to 
atmospheric turbulence 

Let the linear system have the transfer function G(s) which 
relates the output 2(f) (it may be one of many) to the single input 
v(t If x(t) is a transient function, zero for ¢ < 0, then we have 


from Parseval’s theorem? that 


, ’ 1 . 
E [ s*(t)dt = c(w) |\"dw = s(w)t*(w)dw (1 


+ More precisely, the square of the modulus of the amplitude spectrum 


FORUM 825 


where the * denotes the conjugate complex number. £, which 
is easily measured on an anlog computer, may be thought of as 


the ‘‘total energy” of the output, and Z(w) is the Fourier trans- 


o(w J e~'*" 2(t)dl 2 


The output transform Z(w) is obtained from the input trans- 


form 


form by the transfer relation 
Siw X(w)G(tw 3 
hence 
, l : cae 
E X (w) |? |G(tw) |"dw } 
T ( 


When the input is a random function, with power spectral 
density ¢(w), the well-known relation for the mean-square re- 


e «x 
2? 2 o(w) |G(iw) |\*dw 5 
0 


Comparing Eqs. (4) and (5) we see that E is exactly equal to 2? if 


sponse is 


X(w) |? 21o(w (6 
Thus if the function x(t) is such that Eq. (6) is satisfied, then the 
integral of the square of the response to x(t) is exactly equal to 
the desired mean-square response to the given random input 


The Input Functiont ag(t 
A suitable function for the transient gust intensity is 
well = Ae 7 (4 
from which the input function is 


ag(l —(A/uo)e (S 


This gust function has the transform 
dl A 
W,(k) = J a,(lje— "dl A/(y + tk) 
0 


whence 


t For further details of the method of application of this input function 


see Refs. 5and7 
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Fic. 1. Comparison of the spectrum function derived herein 
with the commonly-used function ® (Q) 
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The corresponding power spectral density, from Eq. (6), is 
o'(k 1/2r[A?/(y? + k?)] (10 

The one-dimensional spectrum function commonly used for 
atmospheric turbulence is (Ref. 1, Eq. 2.6 

$(%) = (oh /2r)} (1 + 3(LQ,)2]/[1 + (LQ,)?] 44 (11 
Since Eq. (7) has two unspecified constants, A and y, we can 
impose two conditions on it. We choose, from those which are 
possible, (i) the spectral densities given by Eqs. (10) and (11 
shall be equal in the limit 2; ©; and (ii) the intensity of tur 
bulence (mean-square downwash) given by Eqs. (10) and (11 
From these conditions we get® 


A =6 V3c/2L (12 


and + 3¢/4L (13 


shall be the same 


whence 


ag(l o/twV (3c/2L) exp | —(3c/4L)b) (14 


The power spectral density corresponding to Eq. (10), referred to 


Q2,, 1S 


p'(Q1) o'(k dk /dQ, 
which, after substitution for A, y, and k becomes 


$'(%) = (Be2L/2r)}1/[2.25 + (LA%)*]} (15) 


The two functions ¢ and @’ are compared in Fig. 1. It is seen 
that the difference between them is small except at the lowest 
wave numbers (longest wavelengths). It is precisely here that 
the experimental evidence is most inconclusive (Fig. 10.5 of Ref. 
6) and it appears that there is no good reason to prefer one of these 
functions to the other. It is concluded that $’(Q) gives a good 
representation of atmospheric turbulence, and hence that Eq 
(14) describes a suitable transient. 

The method described above has recently been applied to a 
problem of control optimization’ and found to be very con- 


venient. 
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Comments on 
‘‘A Sublayer Theory for Fluid Injection’’ 


Bernard M. Leadon* 

General Dynamics Scientific Research Laboratory, 
San Diego, Calif. 

March 6, 1961 


SUMMARY 
Three comments on Professor Turcotte’s recent paper! should 
be made: (1) it is not evident that linearity of the velocity pro- 
file in the semilog plot at once justifies Prandtl’s mixing-length 
assumption for the turbulent-injection boundary layer; (2) one 
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Fic. 1. Typical shear-stress and velocity distributions with 
and without air injection into a low-speed, turbulent boundary 
layer. The outer region is affected much as in the case of adverse 
pressure gradients, but the inner regions behave uniquely under 
the influence of the inertia of the injected air 


can achieve a correlation equivalent to that presented without 
assuming anything about the form of the sublayer, buffer-layer 
velocity distribution; and (3) the similarity parameter used is, 
on present indications, inferior to one already in common use.?**7 
These points will be elaborated upon in the following discussion 


DISCUSSION 


ges EQUATION of motion combined with the continuity equa 
tion may be integrated formally to obtain 


r/pU? Tw/pU? + (vy/U)z2 4 
. . 
(0/O0x) | sdy — 2(0/0x) | gdy (1) 
e e 
where s = u/U, U is the free-stream velocity, subscript w refers 
to the wall (y = 0), and the indefinite integrals are on the in- 


terval (0,y). When the upper limit tends to infinity (or 6), the 
magnitude of 7 diminishes to zero, and the equation assumes the 
familiar form of the von Kdrmdn momentum-integral equation 
including the effect of injection or suction. For sufficiently 
large values of x or small values of y the last two terms may be 
neglected. If v» = 0, this approximation indicates that r = ry 
near the wall. 
portion, the ‘overlap’ 


Thus it follows that the existence of a linear 
’ region, in a plot of u versus log y suggests 
that the mixing length, / V 1/p/\0u/Oy, may be taken to 
be simply proportional to y as a useful postulate. This justi 
fication fails, of course, when v,, + 0 

Still 7 may possibly be constant in the overlap region owing to 
contributions from the last two terms of the equation. This 
would be true if the magnitude of the injection velocity were in- 
significant there. To resolve this uncertainty the complete 
equation was integrated numerically® by direct use of the avail- 
able experimental data’ aided by a previous correlation.‘ The 
results are shown in the figure for a typical blowing condition 
and for zero blowing. Also shown are the corresponding velocity 
distributions which exhibit the parallel overlap regions noted 
earlier. One may draw three conclusions from the curves: (1 
the effect of fluid injection upon 7 and uw is not confined to the 
sublayer; (2) 7 is not constant with y when v, > 0; and (3 
Prandtl’s mixing-length assumption should be modified to ac 
count for the injection velocity. (This last conclusion has also 
been acted upon.®) 

Skin-friction data for various injection rates may be correlated 
either by an application of dimensional reasoning or by direct 
use of the above equation without making any assumptions as 


to the nature of the turbulence. Dropping the last two terms, 
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which in the numerical calculations were found to be negligible 
for YV two/p/v < 500, and dividing each remaining term by 7. 
pU?, where subscript 0 refers to zero blowing, we may designate 


by #ao* that value of u* for which 7 = 7, , and finally write 
7/ Tw 1 — Ugo Vu (2 
where ( )y* indicates ( )/V rw/p. This relationship sug- 


gests that v.9* may be an adequate similarity parameter. Asa 
limiting condition, if even incremental injection immediately 
affects the entire sublayer and buffer layer, as the profile data 
suggest, the limiting value of ua)* will be about 13. To the ex- 
tent that this is true we may say that our theory is complete for 
small injection (or suction) rates without need to determine arbi 
trary constants 

Now to fit the equation to the data for all injection rates we 
may try in turn a number of simple assumptions as to the varia- 
tion of uq* with injection rate. For example, we may seek an 
intermediate ‘‘friction velocity’’ which will always be in constant 
proportion to u%@. Among the simplest choices then are: (a) 
Ut = ki, (b) V Two/tw Uaot = ke, and (c) V2rw/(tw + tw) 
“ao? = k3, where ki, ko, and k; are constants. These alternatives, 


when inserted into Eq. (2), yield, respectively: 


T/Tw = 1 — Rivuo* (3) 
Tw/ Tw lL — keV Tw/ Two Vwo* (4) 
ind Tw/ Tw 1 — ksV (rw + tw)/ Two Vwo * (5) 


Iu comparing these possible forms with the experimental data it 
uppears that Eq. (5) with k; = 15 is satisfactory. For all prac 
tical purposes it is equivalent to Turcotte’s Eq. (18). 

However, there is a valid theoretical objection to the use of 
vyo* alone as the injection parameter. This arises in the desire 
to extend such correlations to dissipative boundary layers 
wherein the fluid injection mechanism may be used for heat 
protection. In such cases, if the heating rate is to be constant, 
the blowing rate must be distributed approximately in the same 
way as the skin friction.2-° The proper similarity parameter then 
seems to be Ly)*?.9*. A number of experiments in high-speed 
flows have tended to confirm the usefulness of this blowing param- 
eter,’ but the evidence is admittedly not yet conclusive 

On the other hand, in support of the parameter vy .* one has 
only the evidence of the foregoing theory and the correlation pre- 
sented by Turcotte, plus an argument from the ‘‘law of the wall,”’ 
the one relying entirely upon wall-friction data the accuracy of 
which at high blowing rates may be questioned, and the other 
appealing again to the hypothesis that all injection effects are 
confined to the sublayer. Observing the fact that the momen- 
tum defect region is also strongly affected by mass injection, pre- 
pares one to accept the appearance of the free-stream velocity in 


a satisfactorily chosen injection parameter. 
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Author's Reply 


D. L. Turcotte 

Professor, Graduate Schoo! of Aeronautical Engineering 
Cornell University, Ithaca, N. Y 

April 3, 1961 


M”" AUTHORS who have studied the problem of fluid injec- 
tion into the incompressibie turbulent boundary layer have 


used the momentum equation valid near the wall 
T Tw + plwtt l 


The real question, however, is when and how to cut off the solu- 
tion to Eq. (1 This author proposed! the hypothesis that the 
shear stress in the fully turbulent portion of the boundary layer 
is unaffected by injection. The rationale for this hypothesis is 
that the mechanism of turbulent shear dominates the flow in the 
turbulent boundary layer and that injection of fluid at the wall 
cannot affect the shear stress in the fully turbulent portion of 
the boundary layer. The reduction in shear stress predicted 
from this hypothesis was shown to be in agreement with the 
measurements of Mickley and Davis.? It was also shown that 
this hypothesis is consistent with the measured velocity profiles 

Of course the hypothesis of constant shear stress is subject to 
experimental verification. Leadon* has reflected the velocity 
profiles given in Ref. 2 through the integrated momentum equa 
tion by numerical integration to obtain local values of the shear 
stress. The results are not in agreement with the above hy- 
pothesis. However the question arises whether the computations 
are sensitive to changes of velocity in the transition region where 
experimental points are in doubt and where momentum adjust- 
ments must take place according to the theory given in Ref. 1 
If the numerical computations in Ref. 3 are accepted, then the 
basic assumption made in Ref. 1 is not tenable and an alternative 
formulation of the problem must be found 

If the effect of injection is restricted to the sublayer and turbu 
lent core of the boundary layer, dimensional arguments lead 
directly to the blowing parameter v0 * as previously proposed by 
van Driest.‘ If the injection is not restricted to these regions, 
dimensional arguments do not lead to a particular choice for the 
injection parameter. At the present time experimental results 
are not available over a sufficient range of Reynolds numbers to 
verify the validity of any injection parameter 
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Determination of Final Temperature in a Gun 
Tunnel 


Bo Lemcke 
The Aeronautical Research Institute of Sweden 
January 19, 1961 


Sects NOW it has not been possible to achieve the high tem- 
peratures that theory has predicted for gun tunnels. Fur- 
thermore a large temperature drop during a run has been noticed 
Both of these effects have been more pronounced when the final 
piston velocity has been increased. No satisfactory explanation 
for these phenomena has been presented up to this time 

It is shown in this note that the method of temperature deter- 
mination used previously’? has not included a significant factor 
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which greatly affects the estimate of the final temperatures at- 
tainable. Calculations which include this factor yield values 
which agree realistically with the currently available experimental 
data. 

When the piston is accelerated down the barrel a shock wave 
is formed in front of it whose strength increases with piston 
velocity. The problem of the piston and shock history has been 
thoroughly investigated previously.* 4 Due to the piston mass 
neither the piston nor the shock wave will reach their limiting 
speed (that of the contact surface and the shock wave in a shock 
tube) before the end of the barrel. This means that the shock 
strength, p;/po, will increase through the whole barrel. 

If one studies the heating of a certain air layer, Ax; (Fig. 1), its 
temperature at the point (/;) is due to (1) compression through a 
shock wave to the pressure p;, and (2) an isentropic compression 


to £;;. The temperature 77; is then 


Tr pi Way + D/(y — VI + (i/ po) (2 = 
T pol + (vy + 1)/(y — 1)) (pi/ po) \ pi (1) 


When the shock wave is reflected from the end of the barrel, a 
new temperature jump occurs which is given by 
T1 Pui Wy + 1)/y — 10) + (br /p1) 
Ti; pis V+ My + D/(y — 1) bi / pr, 
where py1;/p1; is essentially constant between the end of the barrel 
and the piston and is given the value py;/pr. 

The gain in temperature through this shock wave is much less 
than in the first one. The same is true to a greater degree for 
each successive reflection of the shock wave. In the calculation 
of the final temperature all but the first two shocks may be 
neglected. From the conditions behind the second shock wave 
an isentropic compression may be assumed to the final pressure 


without introducing appreciable error. Thus 


y-1 
= -| F ] " (Tu/T)M\T:/To) = 
T» (pu/ pr) pi 
ea, yore pu I+y, Pb 
px a) 7 (2) y—-1 pr (2 y—-1 po (a 
pu/ pr pr) rtd pup) rth pi - 
y—-1 pr y—l1 po 


This expression cannot be evaluated without knowledge of the 
shock history in the barrel. Fig. 2 shows the minimum tem- 
perature (p;/po) = 1) and maximum temperature (p;/po = 
| Pi/Polshock tube) With a reasonable connecting curve, for helium 
or air as the driving gas. The maximum point corresponds to 
the end of a very long barrel, or to the use of a piston whose mass 
approaches zero. Inspection of Fig. 2 shows that there always 
appears a region immediately in front of the piston where the 
temperature increases. The length of this region depends 
on the pressure ratio p/p, the piston mass, and the gas com- 
bination used. All gun tunnels in use now have, as far as is 
known, values of these parameters which correspond to a barrel 
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end near the point E in Fig. 2. This means that a large tem- 
perature gradient exists between the piston and the end of the 
barrel, and also that the mean temperature is much lower than 
the maximum possible. 

Earlier calculations of final temperature have been based en- 
tirely upon the strength of the first shock wave immediately 
before it is reflected from the end of the barrel. 

In view of the preceding theoretical considerations the gun 
tunnel must obviously satisfy the following condition if higher 
temperatures are to be obtained. The piston has to be acceler- 
ated to its final velocity in some small fraction of the barrel 
length, say 1/5. Then the shock strength will be approximately 
constant during ~70 percent of the barrel. To fulfill this condi- 
tion several possibilities are immediately evident. The most 
obvious solution is that of reducing the piston mass. This, 
however, is a difficult problem since the mass needed is about 
1 gram for a 40 mm diameter piston and a barrel length of about 
3m. Secondly, an increase in the length of the barrel is possible 
The advantage of this, however, is limited by the growth of the 
boundary layer, and the corresponding shock attenuation 
Finally, an evacuated section for acceleration of the piston could 
be introduced in front of the main barrel. In this case it is 
necessary to have a second diaphragm located in front of the 
piston. This diaphragm, however, will interfere with the air- 
flow unless a suitable system for removing it is found. Using a 
burnable material is perhaps the easiest solution 

An accelerating section will be constructed and tried at the 
FFA to see if a noticeable change in temperature and temperature 


drop during a run will occur. 
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On the Shock Shape and Pressure Distribution 
About Blunt-Nosed Cylinders 
Using Blast-Wave Theory 


R. J. Swigart 

Senior Aerodynamicist, Lockheed Missiles and Space Division, 
Palo Alto, Calif. 

November 17, 1960 


_— PURPOSE of this note is to present various general forms 
of the equations for shock-wave shape and surface-pressure 
distributions on blunted cylinders in hypersonic flow as obtained 


from blast-wave theory. Restrictions on Lees’ and Kubota’s 
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BOW SHOCK 





Fic. 1. Flow about a blunted cylinder in steady flight 


equations! are pointed out, and it is shown how they are special 


forms of the general case. 


SYMBOLS 


Cp, = body nose-drag coefficient 

d diameter of cylindrical afterbody (Fig. 1) 

M « free-stream Mach number 

p surface pressure 

Pm = stagnation-point pressure 

Pa = free-stream pressure 

R = shock radius measured perpendicular to body centerline (Fig. 1) 


distance from body nose measured along the centerline (Fig. 1) 
= ratio of specific heats 

A. Sakurai’s second-order solution to the cylindrical blast- 
wave problem has recently been extended to include third-order 
terms.2 Using this third-order solution, identifying the nose 
drag of a blunted cylinder with the energy per unit length of the 
original line charge in the blast-wave problem, and employing 
the hypersonic unsteady analogy results in the following ex- 
pressions for shock-wave shape and surface-pressure distribution : 


1.5730 <x 


M.24/ Cp, 4 
5.8037 (: 2 :, 
M..*Cp, \d M4) 


P _ 0.087 2/ Co, + 0.405 + 0.150 
p x/d 


0.79514 (Cp,)“* 4] ~ E 
= U.iYoO (Cp,,) ‘\ 
d 


( 


x/d 
M.2V/ Co, 


Eqs. (1) and (2) are more general than the corresponding second- 
order equations reported by Lees and Kubota! in that they de- 
pend parametrically on body nose-drag coefficient. Lees’ and 


Kubota’s equations 


R as 1.62 x 
= ().78 1+ (3) 
d Na | M.2d 
p . M.? 
= ().0665 + 0.405 (4) 
p x/d 


have embodied in their derivation a nose-drag coefficient for a 
= 1.4 as determined by modified 


‘.4 


» ¥ 


hemisphere at M/.. = 


. l ( Pn 1 
c yy = ; + 5 
U a p (9) 


Newtonian theory 


For all Mach numbers above and slightly below 7.7, the neglect 
of nose-drag variation with Mach number is justifiable. How- 
ever, for nose geometries that differ considerably from a hemi- 
sphere, a value of Cp, corresponding to that geometry should be 
used in Eqs. (1) and (2), rather than Ea. (5). 

A further restriction of Eqs. (1) to (4) is the use in their deriva- 
tion of a constant specific-heat ratio, y, of 1.4. Since the solu- 
tion to the higher second- and third-order problems is numerical, 
the y dependency cannot be made explicit. By using Sedov’s® 
analytical solution to the first-order problem (MM. —~ ©), how- 
ever, the parametric y denendency may be obtained. 
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The results are: 


R (: Cp ) ’ Ix 
= | (6) 
d 8 aly) \ d ; 


and 
b mw Cy 
/ Ki Y) (7) 
p \ sSaly) 
where 
2 
. \ a) } ‘i : 
K(y7) “SEER (8) 
(9) 2-7 


and a@ (y) is given graphically on p. 231 of Ref. 3 
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A Conical Static Pressure Probe 


H.R. Vaughn 

Supervisor, Aerodynamic Research Division, Sandia Corporation, 
Albuquerque, N. M 

January 16, 1961 


A CYLINDRICAL PROBE for measuring the static pressure of a 
supersonic stream typically must have a tip-to-orifice dis- 
tance of 30 diameters to accurately measure static pressures 
through the Mach number range of 1.2 to 5.0. This results in a 
long, flexible probe (usually 50 diameters), delicate to handle in 
wind tunnels and quite impractical for flight 

A conical probe only 19 diameters long which measures static 
pressure with good accuracy within the range 1.1 < M < 4.54 
(and probably higher) is shown in Fig. 1. Its strength is such 
that flight use is practical up to the point where aerodynamic 
heating becomes excessive 

The downstream shoulder, say a 30° included-angle cone, 
should be 50 boundary-layer thicknesses downstream of the 
rear orifices, or, for turbulent flow, perhaps only an inch down- 
stream. The conservative orifice-to-shoulder distance of 1.5 
in. has proved satisfactory in both wind-tunnel and flight tests 
and further shortening of this dimension seems possible 

This probe was developed from three-dimensional method of 
characteristics considerations, and extensively tested in the 12 
and 20 in. tunnels at the Jet Propulsion Laboratory, the Naval 
Supersonic Tunnel at the Massachusetts Institute of Technology, 
and in flight. Fig. 2 shows the measured static pressure error as a 
fraction of both the dynamic pressure and the static pressure 
Fig. 3 shows a comparison of wind-tunnel and flight data. From 
the probe wind-tunnel performance data shown in Fig. 2, it may 
be seen that the probe indicates a static pressure slightly higher 
than the true static pressure 

Flight tests have shown that serious errors in the probe reading 
develop below M = 1.1, amounting to several percent of static 
pressure. These errors result from the passage of the body bow 
shock wave over the orifices at transonic Mach numbers as would 
be expected. No wind-tunnel tests have been conducted above 
M = 4.54 
least up to M = 7.0 

Tunnel techniques used in this probe development finally 
yielded testing accuracies in reading probe static pressure of 
approximately 0.02 percent g. This accuracy is substantially 


However, theory indicates no additional errors at 
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better than the writer has seen previously. It was obtained 
by the following procedures 


1. All pressure lines were exhausted for 24 hours preceding 


a test, and tunnel measurements were taken only after the tunnel 
had run 10-20 min and had ‘‘settled down.”’ 
2. The Mach number assumed correct was computed using 


al 
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the wind-tunnel stagnation pressure measured in the settling 
chamber and the pitot pressure measured in the test section at a 
point on the tunnel centerline. 

3. The “correct”? static pressure was computed from the 
above Mach number and the settling-chamber stagnation pres- 
sure using isentropic flow relations 

4. It was found that serious errors arose if the tunnel was shut 
down to install the static probe (in the same location as the pitot 


probe in step 2) and then restarted and reset to supposedly the 
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Fic. 3. Comparison of wind tunnel average reading, and 
flight tests. The flight data include some lag, have known camera 
and telemetry errors, and are at a higher Reynolds number 
than the tunnel work. 
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same stagnation pressure (and Mach number). 


errors were associated with minute thermal distortions of the 


tunnel structure. To avoid these difficulties, a traversing gear 


was arranged which replaced the pitot probe with the static 
probe in a few seconds with the tunnel running, and a static 


pressure reading was then taken 


in intermittent tunnels, it is not believed that their pressure 


could be held constant enough for measuring to the accuracy 
desired in these tests 

5. Pressure measurements were made using precision micro- 
manometers with mercury for a fluid 

6. Experience with other similar probes through an angle-of- 
attack range indicated that the amount of flow angularity in the 
tunnel was of no consequence. In some flight applications, the 
probe is mounted on a swivel to eliminate the angle-of-attack 
error 
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The Qualitative Dynamics of Atmospheric 
Entry and Re-Entry} 


John D. C. Crisp* 


Research Associate Professor of Applied Mechanics, 
Polytechnic Institute of Brooklyn, Brooklyn, N.Y. 


January 26, 1961 
F” PLANE Keplerian motion in a central force field 


2E = (e* — 1)/P = V? — 2/R 


E, VP RVo being the energy and moment of momentum per 
unit mass of a small body and V, R the nondimensional speed 
and radial coordinate. The first relation represents a family of 
hyperbolas in an E-P plane, each member characterized by a 
fixed eccentricity e—see Fig. 1, in which f = 2(E + 1). During 
zero-lift motion through the atmosphere, E and P become time- 
dependent and a characteristic point in the f-P plane moves from 
an initial position such as A defining an approach orbit outside the 
sensible atmosphere. The differential relation! isdE/dP = d(E + 
1)/dP = (E + 1/R)/P independently of the form of the drag force 
Now R is within 3 percent of unity for motion within the atmos- 
+ This note arose out of work supported by the Aeroballistics Division 
Marshall Space Flight Center, Huntsville, Alabama, under Contract No 
DA-30-069-OR D-2639. 
* Now, Senior Lecturer in Engineering, Monash University, Australia 
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phere so that 1/R ~ 1; EF + land Parethuslinearly related. For 
impact in one pass, the impact condition? P > 0, f—~>2(1 —1/R) < 
0 is given by B. The atmospheric trajectory is then represented 
ipproximately by motion, at some undisclosed speed, of the char- 
icteristic point along the straight line AB. A minimum e is 
reached at a point such as C which corresponds closely to circular 
speed, f ~ 1. A steeper initial entry, D, results in a larger mini 
mum e; a shallower entry, E, can lead to a minimum e¢ of zero and 
instantaneous circularization at F. In multiple-pass entry from 
A, the characteristic point descends to a new stationary point 
The final 


pass to impact, H — B, always contains the minimization of « 


G upon first exit and say to H upon the second exit 


Entry from J results in essentially the same behavior as that 
from A apart from the number of passes necessary; the similarity 
parameter is (f/P A minimization of e cannot occur for e > 1 

This behavior is confirmed qualitatively by numerical solu 
tions! of the exact equations of motion. As a quantitative inter- 


pretation it clearly fails. Entry from L, for example, would 


e>i 


2(e+1) 


f 


Fic. 1. 


ENERGY PARAMETER 





MOMENTUM PARAMETER P=R*V, 


require negative eccentricities prior to impact. Again, the line 
slope is insensitive to realistic changes in the initial pair f,P 
The linear approximation corresponds, in fact, to rectilinear 
motion (of the body) in the absence of a gravitational force; and 
furthermore, the drag force function is unspecified. The approxi- 
mation is best in the neighborhood of minimum e for initial en- 
tries not too steep 
Other simple approximations are easy to develop: 


d(E + 1/R)]/dP (E + 1/R)/P™[1 + (sin y)/A] 


yields a linear relation between E + 1/R and P provided y, the 
flight path angle, is small compared with the total atmospheric 


This is true close to peak deceleration; for realis- 
< 0.0020, at the first peak 


acceleration A 
tic hyperbolic entry the ratio sin y/A 


acceleration 


More generally, define the variables f’ V2*R", P’ V2R™ 
where 1, n are undetermined constants: df’/dP’ f'/P' & the 
excess of speed over circular is V?R — 1 = (m — n)(2 — m + 


n)—'. In an f’ — P’ plane the path of the characteristic point 
can be approximated by a series of straight lines dependent on 


the body speed range 
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A Note on the Solution of the Boundary-Layer 
Energy Equation With an Energy-Source Term 


E. Dale Martin* 
NASA, Ames Research Center, Moffett Field, Calif 
February 4, 1961 


IT NUMERICAL SOLUTIONS to the laminar compressible boundary- 
layer energy equation containing an energy-source term, an 
iteration procedure is frequently required in order to solve the 
“two-point boundary-value problem.’’ The numerical integra- 
tion is begun at the wall with the wall temperature or enthalpy 
given as a boundary condition, and with the enthalpy gradient 
at the wall being repeatedly guessed until the proper value is 
found such that the enthalpy at the outer edge of the boundary 
layer approaches a designated external enthalpy. This method 
has been used, for example, by Howe.! 

A simple method is herein described in which the proper wall 
enthalpy gradient can be calculated, so that the condition on the 
enthalpy at the edge of the boundary layer is satisfied, before 
the numerical integration of the energy equation is carried out 
Thus, by use of this method, one eliminates the iteration process, 
which is at least cumbersome, and which often hampers con- 
siderably the integration process 

Howe! presents an exploratory analysis of the possibility and 
value of shielding the stagnation region of a blunt body against 
radiation heat transfer from the high-temperature gas in the 
shock layer by injecting an absorbing gas into the boundary 
layer. This method of radiation-heat shielding by transpiration 
of an absorbing gas is shown to be potentially very effective 

Because of current interest in the problem studied in Ref. 1, 
this note will deal explicitly with implementation of the solution 
in that analysis. The reader will note that a similar procedure 
could be applied to the numerical integration of other forms of 
the differential equation, with other boundary conditions, where 
an iteration procedure might otherwise be required because 
sufficient boundary conditions at the same point are not given 

The reader is referred to Ref. 1 for discussion of the various 
simplifying assumptions and approximations used in that analysis 
and for the partial differential equations describing the laminar 
compressible boundary-layer flow of a binary mixture of gases 
in the presence of a radiation field. The partial differential equa- 
tions can be transformed into a set of ordinary differential equa- 
tions by defining the following dimensionless functions of the 
similarity variable n: 

f’(n) = u/te, Z(n Tie | 


lin) = (L/jePebtelero’ V 2s - Win pi pf 


where the independent variables s and y are related to the physical 
space variables, x (distance along the body) and y (distance from 


the body surface), by the Levy transformation: 


x vy 
u 
$ = Pell eltel’y?"AX; n oa ph | p dy 2 
a 0 av zSC ~/ t 
and where u is the velocity parallel to the surface, 7 is the total 
enthalpy A + (1/2)u?, J is the radiation intensity (integrated 
over the wave length), p is the density, yu is the viscosity, 7p is the 
radius of cross section of a body of revolution, and C pi/ Pelte iS 
respectively, zero and 


assumed constant The exponent 1 is, 


unity for the planar and axisymmetric cases. The subscripts , 
and ,;, respectively, indicate flow exterior to the boundary layer 
and the foreign absorbing gas in the boundary layer. The set 
of governing equations and boundary conditions of the problem, 


as developed by Howe,! then become 


er te 0 (3) 
g” + Pr fg’ Pri (4) 
Ww" Sc fi 0 (5) 
y aWl = 0 (6) 


* Aeronautical Research Scientist 





832 JOURNAL OF THE AEROSPACE 


where, at 7 = 0: 
ie ff am g=g, =0 l 
W = W., " = W', = Sc full — Ww) 


and,asn—~ ©: 

f'>1, gl, Ll—l, (8) 
The constants Pr and Sc are, respectively, the Prandtl number 
and the Schmidt number; a = KV (paeC)/[(n + 1)8] where 
K is an absorption coefficient, 8 is the external stagnation point 
velocity gradient; and W,, is the following function of f,, and Sc, 


derived by Howe:! 


: \ ” \Se/C. ” pee / 
W. = j1- [u w)8°/Se fof, Uf an | 
/ ; 0 \ 


[where f’,, = f”(O) is determined by f,,; see below]. 


The system of equations, (3) through (6), can be simultane- 
ously integrated numerically. However, the boundary condi- 
tions must all be known at a single point, that is, at 7 = 0 if the 
integration is to be accomplished without an iteration. For 
integration of Eq. (3), values of f"(0) = /f"» have been de- 
termined,? and may be specified, so that the condition f’(«) — 
1 is satisfied for various values of fy. In order to find the correct 
1(0) = 1» to use for Eq. (6), Eq. (5) can be integrated simultane- 
ously with Eq. (3) from » = 0 to © to obtain W vs n, and hence 


to obtain ‘ W dn for use in the following expression obtained 
0 


lw/le = exp [-«f' W dn | (10) 
0 


At this point only Eq. (4) remains a ‘two-point boundary- 
value problem.’”’ The purpose here is to illustrate a simple 
method for calculating g’(0) = g’» so that the condition g( ©) — 1] 
is satisfied, and hence to provide both boundary conditions on 
Eq. (4) at » = 0. Use the fact from Eq. (3) that f = —f’’’/f” 
and define @ = @(n) by 


a= (f")—"™, 0’ = Pr fo (11) 


from Eq. (6): 
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Eq. (4) may then be multiplied through by @ and written as 
(0g’)’ = Pr ol’ (12) 
and integrated once, using the definitions g’(0) = g’, and @(0) = 
6,, to obtain 
0g’ — Owg'w = PrG (13) 
where G = G(n) is given by 
Cat) atx = 0,G =0 (14) 


Eq. (13) may then be multiplied by (1/0) and integrated, using 


the condition g(0) = gy =~ O, to obtain 


g = Owe’ wl\(n) + Pr He(n (15) 


g 

where dH,/dn = wor (ats = 0, H, = 0 

and dH2/dy = (f")"'G Git» = 0; ZH, = 0 (17) 

The condition g( ©) — 1 may be applied to Eq. (15); the result is: 
gw’ = [(f"w)"/Ai(@)| [1 — Pr H2(@)] (18) 


which may then be used as a boundary condition at » = 0 to re- 

place the condition g(“) —~ 1. The calculation procedure is as 

follows: Integrate Eq. (3) from 7 0 to © to obtain (rr ya 

0 * 

dn and evaluate W, and W’,. Integrate Eqs. (3) and (5) simul- 
el 

taneously from » = 0 to © to obtain W’ dyn and evaluate 
0 

l, from Eq. (10). Integrate Eqs. (3), (5), (6), (14), (16), and 

(17) simultaneously from » = 0 to to evaluate H;(@) and 

H,( © ) and hence to calculate g,,’ from Eq. (18 

Now Eqs. (3) to (6) are easily handled by simultaneous ma- 

chine integration, without iteration, starting at » = 0 with the 

quantities f., f’. = 0, f’w, Ze = 0, g’«, We, We, and l,, which are 

either specified or calculated as described above 
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Longitudinal Dynamics Continued from pave 78s) 





is the flight path angle (positive for a descending path), 
y=a-0 


The quantities on the right-hand side of Eq. (A.1) are 
those which are excited in the phugoid mode. They 
are shown in the appendix to reference 8 to be related 
to one another in this mode as 


th = e€ COS pnt r= —€ COS Oy». 

vy = esin &nf 
where e¢ is the small eccentricity of the elliptic orbit. 
Thus we have the picture of the pitching motion being 
driven by the trajectory perturbation. When the two 
frequencies are very far apart, the coupling is weak. 
When thev are equal, there is a resonant build-up of the 
pitch angle @. 
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makes it difficult to establish the identity of the author. This is 
particularly important for large annual indexing and abstracting 
services. All titles and degrees or honors are omitted. The 
name of the organization with which the author is associated 
should be placed after his name on a separate line. The date on 
which the paper is received will be inserted by the Editor. The 
author’s title or position should be indicated in a footnote. 

SUMMARIES OR ABSTRACTS: An abstract to be printed at the 
beginning should accompany each article. It should be adequate 
as an index and asa summary. It should contain a statement of 
major conclusions reached, since summaries in many cases con- 
stitute the only source of information used in compiling scien- 
tific reference indexes. Abstracts printed in other journals, espe- 
cially foreign, in most cases consist of summaries from printed 
papers. The summary should explain as adequately as possible 
the major conclusions to a nonspecialist in the subject and should 
contain from 100 to 300 words, depending on the length of the 
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at frequent intervals. The work of editorial preparation will be 
simplified by the author’s provision of many subheadings. 
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little technical interest and not showing advances in general 
practice. Too detailed tabular matter (general results of such 
tables may be included in the text). Lengthy descriptions of 
materials or processes or of preliminary experiments or theories 
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ILLUSTRATIONS: Illustrations should accompany manuscripts, 
and each should always be referred to in the text by number. 
Drawings or graphs should not be larger than 12 X 16 inches, 
and must be made with jet black India ink on white paper or 
tracing cloth, the latter being preferred. Do not use typewriter 
The smallest lettering on 8 X 10 inch figures should 
be no less than '/, inch high. Cross-section paper (white with 
black lines) may be used, but it should not have more than 4 lines 
per inch. If finer ruled paper is used, the major division lines 
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In the case of finely ruled paper, only blue-lined paper can be 
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be on glossy white paper. Avoid round or oval photographs. 


CAPTIONS AND LEGENDS: Legends or captions must accompany 
each drawing or photograph submitted. If written on the draw- 
ing or photograph they should be placed below and well outside 
the part to be reproduced. Each table should have a caption 
such as Table 1, Table 2, Table 3, etc. Captions should be com- 
plete in themselves so as to make the data intelligible to the 
reader without reference to the text. A duplicate list of captions 
for figures should be included as the last page of the manuscript. 
Use “‘Fig. 1’’ (not Figure 1), “‘Figs. 3 and 4,” etc., in both the text 
and the numbering of illustrations. In the text, “Eq. (1)” or 
“Eqs. (1) and (2)” is used, not “Equation (1).” In captions, 
legends, and in table headings, write all words in full; do not 
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MATHEMATICAL WoRK: Formulas may be typewritten or 
carefully written in pen and ink, the writing to be large enough so 
that it can be marked for the printer. Considerable space for 
marking should be allowed above and below all equations. All 
complicated equations should be repeated on separate sheets with 
adequate space left for marking. The solidus should be used for 
simple fractions appearing within the text. Make all expressions 
clear to the typesetter. Greek letters used in formulas should 
be clearly designated by name in the margin of the manuscript. 
The difference between capital and lower-case letters should be 
clearly distinguished and care taken to avoid confusion between 
zero (0) and the letter (0), between the numerical (one) and the 
letter (ell) and the prime (’), between alpha and a, kappa and 
k, u and mu, v and nu,n and eta. All subscripts and exponents 
should be clearly distinguished. Avoid complicated exponents and 
subscripts. Dots and bars over letters or mathematical expres- 
sions should also be avoided. When it is necessary to repeat a 
complicated expression, it should be represented by some con- 
venient symbol 


SYMBOLS AND ABBREVIATIONS: The symbols recommended in 
the American Standard Association “‘Letter Symbols for Aero- 
nautical Sciences,’”” ASA Y10.7—1954, should be used wherever 
practicable. All symbols should be clearly written and carefully 
checked. Standard abbreviations should be used, and it should 
be noted that most abbreviations are lower case. 
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